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ABSTRACT 

In this study we present the results from reahstic A^-body modeUing of massive star 
clusters in the Magellanic Clouds. We have computed eight simulations with N ~ 10^ 
particles; six of these were evolved for at least a Hubble time. The aim of this mod- 
elling is to examine in detail the possibility of large-scale core expansion in massive 
star clusters, and search for a viable dynamical origin for the radius-age trend ob- 
served for such objects in the Magellanic Clouds. We identify two physical processes 
which can lead to significant and prolonged cluster core expansion - mass-loss due to 
rapid stellar evolution in a primordially mass segregated cluster, and heating due to a 
retained population of stellar-mass black holes, formed in the supernova explosions of 
the most massive cluster stars. These two processes operate over different time-scales 
and during different periods of a cluster's life. The former occurs only at early times 
and cannot drive core expansion for longer than a few hundred Myr, while the latter 
typically does not begin until several hundred Myr have passed, but can result in core 
expansion lasting for many Gyr. We investigate the behaviour of each of these expan- 
sion mechanisms under different circumstances - in clusters with varying degrees of 
primordial mass segregation, and in clusters with varying black hole retention frac- 
tions. In combination, the two processes can lead to a wide variety of evolutionary 
paths on the radius-age plane, which fully cover the observed cluster distribution and 
hence define a dynamical origin for the radius-age trend in the Magellanic Clouds. We 
discuss in some detail the implications of core expansion for various aspects of globular 
cluster research, as well as the possibility of observationally inferring the presence of 
a significant population of stellar-mass black holes in a cluster. 

Key words: stellar dynamics - globular clusters: general - methods: iV-body simu- 
lations - Magellanic Clouds. 



1 INTRODUCTION 

As relatively simple objects which are integral to the study 
of many fundamental astronomical processes, massive star 
clusters are central to a wide variety of astrophysics over 
all scales - from star formation and stellar and binary evo- 
lution, through stellar exotica and variable stars, and the 
dynamics of self-gravitating systems, to galaxy formation 
and evolution, with implications for cosmology. In the con- 
text of this wider astrophysics however, it is clearly essential 
that we understand the clusters themselves: how internal 
physical processes in clusters shape their overall character- 
istics (and vice versa), and how individual clusters interact 
with and are influenced by their local environments. Only 
when armed with this knowledge is it possible to disentangle 
cluster evolutionary processes from the specific astrophysics 
under investigation. 



From an observational perspective, we are provided 
with only a limited set of massive stellar clusters which are 
close enough to us that they may be fully resolved using 
presently available facilities (and hence thoroughly studied 
on a star-by-star basis) . The Galactic globular clusters, while 
constituting the closest ensemble, are not ideal for study- 
ing massive star cluster evolution, primarily because they 
are almost ex clusively ancient objec t s with ages ~ 10 — 13 
Gyr (see e.g., [R osenber g et al.|[l99 p: 'Salaris fc Weissll20oi 
iKrauss fc Chab ovcr 20031: iDe Ange h ot al. 2005). Therefore, 
while we are able to precisely measure the end-points of mas- 
sive star cluster evolution, the long-term development which 
brought them to these observed states must be almost com- 
pletely inferred. Fortunately, it is relatively straightforward 
to turn our attention to the Large and Small Magellanic 
Clouds (LMC/SMC), which both possess extensive systems 
of star clusters with masses comparable to the Galactic glob- 
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ulars, but crucially of all ages: 10® i$ t ^ 10^° yr. These 
two nearby galaxies are hence of fundamental importance to 
studies of star cluster evolution, because they are the only 
systems in which we can directly observe snapshots of clus- 
ter development over the last Hubble time using a sample 
of fully resolved objects. 

Some of the earliest studies to take advantage of this 
situation and investigate the structural evolution of mas- 
sive stellar cl usters were those of El s on an d collaborators. 
In particular, lElson. Fall fc FreemanI ijigSTh constructed ra- 
dial brightness and den sity pro files for 10 young cluster s in 
the LMC, w hile IeIsou. Freeman fc Laueil (1989) and Elsoiil 
l|l99ll . ll992l ') extended this study to a larger sample of LMC 
clusters including much older objects. They discovered a 
striking relationship between cluster core size and age - 
specifically, that the observed spread in core radius is a 
strongly increasing function of age. The youngest clusters 
in their sample possessed compact cores with rc ~ 1 — 2 pc, 
while the oldest clusters exhibited a range < rc < 6 pc (cf. 
Fig. [T|. Here, cluster core size is parametrised by the obser- 
vational core radius, Vc, defined as the projected radius at 
which the surface density/brightness has decreased to half 
its central value. 

The advent of the Hubble Space Telescope (HST) 
has allowed this discovery to be re-addressed observation- 
ally in significantly more detail than was possible with 
ground-based facilities. HST imaging can resolve Magel- 
lanic Cloud star clusters even in their inner cores, so that 
star counts may be conducted to very small projected radii 
and accurate surfac e density /brightness profiles constructed. 
iMackev fc Gilmord (|2003al ) obtained structural measure- 
ments from a homogeneous compilation of archival Wide 
Field Planetary Camera 2 (WFPC2) imaging of 53 mas- 
sive LMC clusters spanning the full ago range . We found 
essentially the same relationship as Elson ct alj (|l989l ) - the 
youngest massive LMC clusters possess compact cores of 
typical radius ~ 1 — 2 pc, but with increasing age the spread 
in core radius increases such that the oldes t cluste rs span 
the range < rc < 8 pc. IMackev fc Gilmord (|2Q03lJ) subse- 
quently extended these HST measurements to 10 SMC clus- 
ters, demonstrating for the first time that a radius-age trend 
indistinguishable from that observed in the LMC exists for 
this star cluster system. 

Following these two studies, we were granted HST time 
to conduct a snapshot survey of additional massive LMC and 
SMC star clusters using the Advanced Camera for Surveys 
(ACS; HST program #9891), with the aim of improving the 
sampling of the radius-age parameter space. In all, 31 extra 
LMC and 13 extra SMC clusters were successfully imaged, 
significantly enlarging the sample. Final structural measure- 
ments for these objects are yet to be published (Mackey 
et al. 2008, in prep.); however preliminary results for the 
core sizes are plo tted in Figure U al ong with the WFPC2 
measurements of IMackev fc Gilmord (|2003al lbl) . In obtain- 
ing these new parameters, photometric measurements were 
made using the pipeline described by IMackev fc Gilmord 
l|2004l and IMackev. Pavne fc Gilniord (|2006i'). while radial 
brightness profiles were constructed follo wing procedures es- 
sentiall y identical to those described bv IMackev fc Gilmord 
l|2003al ) but adapted for ACS Wide Field Channel (WFC) 
imaging. Figure [T] also includes the r ecent HST/ACS mea - 
surements of the SMC cluster BS90 bv lRochau et all (|20o3). 
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Figure 1. Core radius versus age for massive stellar clusters in tfie 
Large and Small Magellanic Clouds. This Fig ure includes all clus- 
ters fro m the HST/WFPC2 measurements of lMackev &: Gilmord 
l i2003al lbl) as well as the new preliminary HST/ ACS measurements 
of Mackey et al. (2008, in prep. ) from Program #98 91, and the 
recent measurements of BS90 bv lRochau et all ([SooJ). Core radii 
for several o f the old est compact clusters are upper limits, as in- 
dicated (see IMackev & GilmorQi2003ai ). 



Note that the core radii for several of the oldest, most com- 
pact clusters in Fig. [T] are upper limits, as indicated. This is 
due to severe crowding in the HST imaging, and the possibil- 
ity t hat several of these cluste rs are core-collapsed objects 
(see IMackev fc Gilmord l2003al ). All affected clusters have 
measured rc < 1 pc. 

Figure [1] represents the most complete and up-to-date 
information presently available regarding the radius-age 
trend in the LMC and SMC star cluster systems. The upper 
envelope is very well defined for all ages up to a few Gyr. At 
older times than this, the full range of core radii observed 
in massive star clusters is allowed. In fact the situation is 
even more dramatic than was appreciated by earlier studies. 
Several of the oldest clusters in the new ACS sample fall off 
the top of the diagram: the Reticulum cluster in the LMC, 
with age r ~ 12—13 Gyr and Tc ~ 14.8 pc; and Lindsay 
1 and 113 in the SMC, with r ~ 9 Gyr and ~ 16.4 pc, 
and r ~ 5 Gyr and rc ~ 11 pc, respectively. Hence the size 
range observed for the oldest clusters is < rc < 17 pc. 

These recent measurements of very extended objects are 
consistent with those for several old globular clusters in the 
Fornax and Sagittarius dwarf galaxies - Fornax c luster 1, 
and Terzan 8 and Arp 2 (|Mackev fc Gilmordl2003cr i - which 
also have very large core radii. A number of Galactic glob- 
ular clusters are also known to possess extended cores, as 
seen in Fig. [2] (upper panel) , which shows the core radius 
distribution of the oldest LMC and SMC clusters from Fig. 
[1] compared with that for the Galactic globular cluster sys- 
tem. The observed ranges in rc match well, as do the general 
shapes of the distributions. The main difference is that the 
distribution for Galactic globulars is more sharply peaked 
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Figure 2. Core radius distributions for the oldest (t > 7 Gyr) 
Magellanic Cloud clusters from Fig. [T] (dashed lines) and for 
Galactic globular clusters (solid lines). Upper panel: The full 
sample of Galactic globular clusters with suitable measurements 
of Tc is plotted; objects which are members or ex-mem bers of 
the Sa gittarius dwarf galaxy are excluded. Data are from iHarrij 
lll996l') (2003 on line update) with new measurements as described 
in lMackev fcG ilmorc (20ol) and|ffilker (2006). The distributions 
are very similar in range and overall shape; however the Galactic 
globulars have a sharp peak at small core radii. Note that all Mag- 
ellanic Cloud clusters with core radius measurements which are 
upper limits already fall in the smallest rc bin. Lower panel: If 
only the Galactic "young-halo" subsystem is considered (objects 
which preferentially lie at Galactocentric radii beyond ~ 15 kpc) 
a very much closer match is observed. 



at small Vc- This is not surprising given that a large frac- 
tion of the Galactic globulars reside in the inner Galaxy, 
where tidal forces are expected to rapidly destroy l oosely 



where tidal forces are expected to rapidly destroy l oosely 
bound cl usters. Indeed, following Mackev fc Gilmora (|2004l ) 
(see also lZinn|[l993 : lMackev fc van den BerghlbOOSL if only 
members of the Galactic globular cluster "young-halo" sub- 
system are considered (most of which are located at Galac- 
tocentric radii larger than ~ 15 kpc), the core radius distri- 
bution is an excellent match to that observed for the oldest 
LMC and SMC objects (Fig.[2l lower panel). The end-points 
of structural evolution observed for the Galactic globulars 
appear quite consistent with the end-points of the radius-age 
trend observed in the Magellanic Clouds. 

The simplest interpretation of the radius-age trend is 
that it represents the progression of cluster structural prop- 
erties with tim43- In this scenario. Figures [1] and [2] pro- 
vide striking evidence that our understanding of massive 
star cluster evolution is incomplete, since standard models 
do not predict an order-of-magnitude expansion of the core 



An additional possibility is that cluster formation conditions 
have changed significantly over the past ~ 10 Gyr. However, 
presently available constraints on this proposition are limited, 
and we do not discuss it further here. 



radius over a Hubble time (e.g., iMevlan fc Heggi3 Il997l ). 
Identifying the origin of the radius-age trend is therefore of 
considerable importance for star cluster astrophysics, and all 
related fields in which star clusters play a prominent role. 

lElson et all (|l989l ) discussed the possibility that the in- 
creasing spread in radius with age could refiect inter-cluster 
variations in the slope of the stellar initial mass function 
(IMF). Clusters with flat IMFs possess comparatively more 
massive stars than those with steep IMFs. Consequently, 
they suffer more severe mass loss due to stellar evolution at 
early times, result ing in increased relative expansion. How- 
ever, lElson et al.l ([1989) found that to induce expansion 
along the upper envelope of the observed trend would re- 
quire a very flat IMF slope, and the resulting early mass loss 
would be severe enough to disrupt the cluster within only a 
few tens of Myr. An additional problem with the IMF hy- 
pothesis concerns the time-scale - the severe mass loss phase 
lasts for roughly only the first ~ 100 Myr of a cluster's evo- 
lution. Therefore, it cannot drive significant expansion over 
the full range of ages observed for Magellanic Cloud clusters. 
There is also an increasing body of observational evidence 
that the IMF in young star clusters is more-o r-less invariant 
(see e.g..lKroupall200li;lde Griis et aLll2002cl) . 

IWilkinson etal\ (|2003l ) used A'^-body simulations of 
small star clusters to investigate whether the radius-age 
trend could refiect core expansion induced by populations 
of binary stars, or by time-varying tidal fields such as those 
which clusters on highly elliptical orbits might experience. 
They observed similar core radius evolution for model clus- 
ters on both circular and elliptical orbits, and therefore con- 
cluded that the tidal fields of the Magellanic Clouds have not 
yet significantly influenced the evolution of the intermediate- 
age clusters in these systems. Furthermore, while they found 
that the presence of large numbers of hard primordial bina- 
ries in their small clusters did lead to a degree of core radius 
expansion, the magnitude of the effect was insufficient to 
explain the observed radius-age trend. 

iHunter et al] (|2003l ) suggested that rather than repre- 
senting the results of dynamical evolution, the radius-age 
trend might instead have its origins in a size-of-sample ef- 
fect. They measured a very large sample of Magellanic Cloud 
clusters with masses 10 < Md < 10® Mq and found the sig- 
nature of such an effect in their data. On a log-abscissa plot 
such as Fig. [T] older ages correspond to larger time intervals 
and hence to more clusters forming in each log-time inter- 
val. Since the star cluster mass function decreases steeply 
with increasing cluster mass, this results in the maximum 
observed cluster mass in each log-tim e interval inc r easing 
with age. The clusters in the sample of iHunter et al.l (|2003l ) 
also showed a weak dependence of size on total mass in that 
more massive clusters have larger characteristic radii. Com- 
bined with the size-of-sample effect, this leads to a size-age 
distribution with an upper envelope not dissimilar to that 
evident in Fig.[T] However, it is not clear how applicable this 
argument is to the cluster sample considered in Fig. [T] be- 
cause all these clusters have masses M > 10" M©, and show 
no coherent link between total mass and core radius (e.g., 
Mackev fc G Umore 2003ai ). Indeed, restricting the sample of 



Hunter et all (|2003i ) to clusters with M > 10* M©, their 



relationship between size and mass is no longer evident. 
Hence, a size-of-sample effect apparently does not explain 
the radius-age trend visible in Fig. [1] 
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Finally, iMerritt et alj (|2004l ) examined the formation of 
cores in primordially cusped clusters (i.e., objects which ini- 
tially have Ta ~ 0) due to the presence of populations of 
massive stellar remnants. They used analytic calculations in 
combination with simplified A'^-body models (composed of 
equal-mass non-evolving particles) to show that the orbits 
of the remnants decay due to dynamical friction so that they 
sink to the cluster centres, heating the stellar background in 
the process and turning the cusp into a core. The authors 
also note that further heating of the core may continue over 
a longer time-scale, due to subsequent evolution of the sub- 
system of m assive remnants. Th e rates of core growth de- 
termined bv lMerritt et all l|2004 ) are moderately successful 
in reproducing the observed radius-age trend; however their 
models seem to require a range of initial densities which 
is significantly larger than that found for young clusters in 
the Magellanic Clouds. It is also not clear how their results 
would respond to the introduction of a mass spectrum and 
stellar evolution into the simulations, or the introduction of 
more realistic initial conditions including the possibility of 
primordial mass segregation. 

As demonstrated above, the radius-age trend is indis- 
tinguishable in the LMC and SMC, and the end-points of 
the trend are consistent with the core radius distributions 
of the Galactic globular clusters as well as of globular clus- 
ters belonging to the Forn ax and Sagittarius dwarf galaxies 
jMackev fc Gilmorell2003d ). These galaxies cover a very wide 
range of masses and morphological types, and hence pos- 
sess very different tidal fields and possible external torques. 
This strongly suggests that the radius-age trend is primar- 
ily driven by i nternal cluster process es, rather than external 
infiuences (cf. IWilkinson et al.|[2003l ). To this end, we have 
conducted a series of large-scale, realistic Ai'-body simula- 
tions of Magellanic Cloud clusters with the aim of investi- 
gating an internal dynamical origin for the radius-age trend. 
More specifically, we have examined the infiuence of stellar- 
mass black holes (BHs) , formed in the supernova explosions 
of the most massive cluster stars, on the long-term evolution 
of massive stellar clusters. We have also investigated the role 
played by primordial mass segregation in shaping the early 
evolution of massive stellar clusters. The basic results from 
several of o ur key simulations h ave been outlined in a re- 
cent Letter (|Mackev et al.l [20071 ): in the present paper, we 
describe in detail the complete results of our modelling. 



2 NUMERICAL SETUP 

2.1 TV-body code and initial conditions 

We use direct, realistic A^'-body modelling in order to investi- 
gate the structural and dynamical evolution of massive star 
clusters in the Magellanic Clouds. Simulations of this type 
are a powerful tool for such work because they incorporate 
all the relevant physical processes with a minimum of simpli- 
fying assumptions. Recent technological developments mean 
that it is now feasible to run models with A'^ sufficiently large 
so as to be directly comparable to observed clusters. This 
has a number of advantages, discussed below. 

F or the presen t stud y, we have used the nbody4 
code (|Aarseth|[l999l . |2003| ) in combi nation with a 32-c hip 
GRAPE-6 special purpose computer (|Makino et al.ll2003l ) at 



the Institute of Astronomy, Cambridge. Th is code uses the 
fourth-order Hermite scheme (|Makindll99ll ) and fast evalua- 
tion of the force and its first time derivative by the GRAPE- 
6 to integrate the equations of motion. Close encounters be- 
tween stars, including stable binary systems and hierarchies, 
are integrated via st ate-of-the-art two-body an d chain reg- 
ularization schemes l|Mikkola fc AarsethI 1 19931 . |l998) . Also 
included in nbody4 are routines for modelling the stel- 
lar evolution of both single and binary stars. For single 
stars the se take the form of the a nalytical formulae de- 
rived by iHurlev. Pols fc ToutI (I2OO0I ') from detailed stellar 
evolution models, following stars from the zero-age main se- 
quence through to remnant phases (such as white-dwarfs, 
neutron stars and black holes) . Binary star evolution is cal- 
culated in a similar manne r, following the prescription of 
iHurlev. Tout fc Polj (|2002l ) and allowing for such phases 
as the tidal circularization of orbits, mass transfer, and 
common-envelope evolution. The stellar and binary evolu- 
tion is calculated in time with the dynamical integration so 
that interact ion between the two is sim ulated in a consistent 
fashion ('e.g.. lHurlev et al.ll200lll2005h . The stellar evolution 
routines allow a spread in stellar masses covering the range 
0.1 — 100 M0, so that one can construct any desired IMF for 
a model cluster. In addition, a uniform metallicity for the 
cluster may be selected in the range Z — 0.0001 — 0.03. A 
mass-loss prescription is included such that evolving stars 
lose gas through winds and supernova explosions. This gas 
is instantaneously removed from the cluster. Such mass-loss 
can rapidly alter the gravitational potential of a star cluster, 
strongly affecting its early structural and dynamical evolu- 
tion. 

When constructing the initial conditions for our sim- 
ulated clusters, we were careful to develop models as sim- 
ilar as possible to the youngest massive clusters observed 
in the LMC and SMQj. In general, young massive Magel- 
lanic Cloud clusters possess radial surface brightness pro- 
files best described by three- parameter models of the form 
jElson. Fall fc Freemanll 19871 : EFF models hereafter): 

M'-p) = MO (^1 + ' , (1) 

where Vp is the projected radius, fio is the central surface 
brightness, 7 determines the power-law slope of the fall-off 
in surface brightness at large radii, and a is the scale length. 
These models ar e a subset of the more general family of mod- 
els presented bv lZhad | 1996h . The scale length, a, is related 
to the observational core-radius, Vc, defined here as the pro- 
jected radius at which the surface brightness has dropped to 
half fio, by: 

rc = a(2^ - 1)5 . (2) 

Some of the global properties observed for young massive 
Magellanic Cloud clusters are summarized in Fig. [3] In this 
plot, we have taken th e struc tural parameters Vc and 7 
from lMackev fc Gilmord (|2003al lbl). who constructed surface 
brightness profiles from HST photometry and fit EFF mod- 
els as defined above. We have also taken the central density 

^ Although we again note the possibility that the initial condi- 
tions for massive clusters which formed at high redshift may be 
different to those for clusters forming today. 
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and t otal mass estimates of iMcLaughlin fc van der Marell 
(|2005l ). which wer e computed in a more robus t manner than 
those provided bv lMackev fc Gilmor^ (|2003al lbh. 

All young LMC and SMC clusters are observed to have 
cored (rather than cusped) profiles - even the ultra-compact 
cluster R136 exhibits a small co re (see e.g., the discus- 
sion in iMackev fc Gilmord l2003al . and references therein) . 
Their profil e s are well fit by EFF models with 7 ~ 2.5: 
lElson et al.1 (|l987l l found a median value of 7 = 2.6 and 
a range 2.2 < 7 < 3.2 for their ten young LMC clusters, 
while the larger sample plotted in Fig. [3] covers the range 
2.05 ^ 7 55 3.79 and has a median value 7 = 2.67. Ex- 
cluding R136, the young LMC and SMC clusters typically 
have central densities in the range 1.6 < log po ^ 2.8 (where 
po is in units of Mqpc"'^), and total masses in the range 
4 < logMtot < 5 (where Mtot is in units of M©). R136 is 
the youngest cluster in the sample, with r ~ 3 Myr, and 
also has the greatest central density: log po ~ 4.8. 

Our model clusters are generated such that they ini- 
tially have structural parameters in projection which are 
consistent with those observed for the youngest LMC and 
SMC clusters. This is achieved by selecting stellar positions 
randomly from the density distribution of an EFF model 
with 7 = 3. Each star is assigned a velocity drawn from a 
Maxwellian distribution, where the velocity dispersion a is 
calculated using the Jeans equations assuming an isotropic 
velocity distribution. In applying this generation algorithm 
it is important to know that for the EFF family of models, 
the deprojected density profile is given by: 



p(r) 



po 



(3) 



where po is the central volume density. From this, we can 
derive expressions for the isotropic velocity dispersion as a 
function of radius. The 7 = 3 case is the closest value of 7 
to the median 7 = 2.67 observed for young LMC clusters for 
which the expression for a is analytic (see Appendix IX)) . 

We assign the star s in each cluster a range of masses 
according to the IMF of lKroupal (j2001), which is a multiple- 
part power- law ^{m) oc m~"' , where ^{m)dm is the number 
of single stars falling in the mass interval m to m -I- dm. 
iKroupal (|200ll ') derived his IMF from a large compilation of 
measurements of young stellar clusters, including many in 
the LMC. This is in contr ast with many other wide l y used 
IMFs - for example, the iKroupa. Tout fc Gilmori (119931 ) 
IMF, which was derived from observations of Galactic field 
stars in the local neighbourh ood and towar ds the Galactic 
poles. Therefore, we prefer the lKroupal (|200ll ) IMF for direct 
modelling of Magellanic Cloud clusters. 

We impose a stellar mass range of 0.1 — 100 M0 for 
our Ai'-body clusters. The extremes of this range are set 
by the lowest and highest mass stars for which reliable 
stellar evolution routines are incorporated in nbody4. Al- 
though stars mo re massive than 100 Mp, do form in large 
star clusters (e.g., IWeidner fc Kroupa|[2006l ). this upper limit 
is p erfectly a cceptable for our present models. For exam- 
ple, iMassev fc Hunter (1998) found only ~ 10-20 stars 
with M > 80 M0 in the extreme LMC cluster R136 (de- 
pending on the adopted stellar evolution m odels), while 
the revised calculations in lMassev et al. I (|2004 l2005). which 
incorporate improved spectroscopy and modelling, suggest 
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Figure 3. Properties of the young massive clusters ob- 
served intheLMCand SMC . Structural data are taken 
from iMackev fc Gilmord ll2003al lbl) . while the central den- 
sity (po) and total mass (Mtot) estimates are taken from 
IMcLaughlin fc van der Marell ||2005| ). as discussed in the text. 



significant reductions in these estimated masses. In prac- 
tice, we expect that increasing our upper mass limit to the 
inferred fundament a l max imum stellar mass ~ 150 Mq of 
IWeidner fc Kroupal (|2004l ) would have essentially no dis- 
cernible effect on the global evolution we observe for our 
models. We note that our lower mass limit means that in 
practi ce we only utilise the exponents ai — as in the lKroupal 
(|200ll ) IMF. 

Our adopted IMF and stellar mass range, along with 
the requirement that our model clusters have masses com- 
mensurate with those observed for young Magellanic Cloud 
clusters (Fig. allows the total number of stars in each 
given model to be assigned. Choosing A'^ « 10^ stars results 
in initial total cluster masses of log A/tot ^ 4.75. 

It is only relatively recently, with the advent of special- 
purpose hardware, that it has been possible to follow models 
with such large N over more than a Hubble time of evolu- 
tion. There are several advantages to running simulations of 
this size. First, the model star clusters are directly compa- 
rable in terms of total mass and central density (see below) 
to the massive clusters observed in the LMC and SMC. We 
are therefore now moving into the regime where many of 
the scaling-with-A issues which it has been necessary to 
account for in previous studies when applying the results of 
A-body simulations to t he evolution of real clusters (see e.g., 
lAarseth fc Heggielll99^ ) are circumvented. In addition, with 
such large N, fluctuations in the global evolution of the A- 
body model are reduced to the point where they are not sig- 
nificant. For small- A models, it has been standard practice 
to average the results of a number of simulations to reduce 
such fluctuations, the amplitudes of which increase with 
decreasing A (e.g.. iGiersz fc Heggi3ll994l : IWilkinson et al.l 
I2OO3I : iHeggie. Trenti fc Hutll2006l ) For large- A models, this 
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process is not necessary (see e.g., iHurlev et aLll2005l i^. Fi- 
nally, with A'^ ~ 10^ we are able to perform detailed sim- 
ulated observations of our models. This allows us to derive 
quantities from the simulations which are directly compa- 
rable to the genuine observations of LMC and SMC star 
clusters. As we discuss more fully in Section |31 this step is a 
vital ingredient in the analysis of realistic iV-body models. 

Star clusters in the LMC are observed at galactocen- 
tric radii between ~ — 14 kpc. We therefore evolve our 
model clusters in a weak external tidal field, rather than 
in isolation. This external field is incorporated by impos- 
ing the gravitational potential of a point-mass LMC with 
Mg = 9 X 10^ M© , and placing the cluster s on circular or- 
bits o f galactocentric radius i?g = 6 kpc. I Wilkinson et al.l 
l|2003l ) give a more detailed description of the implementa- 
tion of the external field within nbody4, which is done by 
integrating the equations of motion in an accelerating but 
non-rotating reference frame, centred on the cluster's centre- 
of-mass. Adopting a point-mass L MC is a significant over - 
simplification; however, as noted bv I Wilkinson et al] (|2003l ). 
the gradient of this potent ial is within a factor of two of that 
in the LMC mass model of lvan der Marel et al.l l|2002t ) at our 
orbital radius. More importantly however, our aim is not to 
examine the e ffect of tidal field s on t he evolution of star 
cluster cores - IWilkinson et al. demonstrated that 

external fields comparable to those experienced by Magel- 
lanic Cloud clusters do not result in strong core evolution. 
Rather, we impose an external tidal field so that the gradual 
evaporation of stars from the cluster may be simulated in a 
self-consistent fashion, and the rates of evaporation between 
different models with the same external potential and escape 
criterion may be easily compared. 

In an external potential, the tidal radius of a star cluster 
on a circula r orbit may be estimated from the relationship 
l|Kinelll962h : 



Mci 
3A/g 



(4) 



where Md is the cluster mass. In nbody4 stars are deemed 
to have escaped the cluster when they reach a rad ius 2rf. 
This is a legitimate approximation - for example iHeggid 
shows that although cluster stars may on occasion 
possess orbits which allow them to move far beyond rt and 
yet return to the cluster, in practice the vast majority of 
stars which move beyond a few rt are permanently lost. 
In our models rt is a non-static quantity (since the clus- 
ter mass is monotonically decreasing with time); therefore, 
the instantaneous value is used when assessing the above 
escape criterion. We caution that other different recipes for 
the implementation of tidal fields exist, which can lead to 
signi ficantly different escape rates and thus cluster lifetimes 
(e.g., lBaumgardtll200ll : iTrenti. Heggie fc Hutll2007l ). 

Within the A'^-body code the equations of motion are 
integrated in scaled units such that G = 1 and the virial 
radius and total mass of the cluster are also set to unity 



^ We caution, however, that small-number statistics may still be 
subject to significant fluctuations between simulations - an ex- 
ample in the present work are the properties of ejected binaries, 
as discussed in Sections 14. II and 14. 21 



(|Heggie &: Mathieulll986l ). For a star cluster in virial equilib- 
rium the initial energy in these units is —1/4 and the cross- 
ing time is 2\/2. Given the total mass of the cluster in so- 
lar masses and an appropriately chosen length-scale (which 
determines the conversion from A^-body units to physical 
units) it is simple to obtain the conversion factors for time 
and velocity from A'^-body units to Myr and kms~^, respec- 
tively. 

Since the chosen length-scale sets the conversion from 
A'-body units to physical units, it controls the physical den- 
sity of the cluster and hence the physical time-scale on which 
internal dynamical processes occur. We assume that a model 
cluster initially just fills its tidal radius. The value of this 
radius at time r = 0, determined via Eq.|4]with Md = Mtot, 
therefore sets the ratio between the length scale in A'-body 
units and in physical units (pc). EFF profiles formally have 
no outer bound, so when randomly generating the initial 
stellar positions we only accept stars lying within the esti- 
mated tidal radius of the cluster under consideration. 

The above process determined a length scale of 8.26 pc 
for our model clusters. This results in an initial central mass 
density of logpo = 2.31 and a core radius rc ~ 1.90 pc for 
these objects, values which match well those observed for 
many young Magellanic Cloud clusters (Fig. We note 
that the clusters described here are not in any way mass 
segregated; however we also ran simulations of clusters in- 
corporating various degrees of primordial mass segregation, 
the details of which are described below in Section[2]2] Those 
objects have the same length scale as the clusters described 
here, but smaller core radii and much higher central densi- 
ties, more in line with those of the very young LMC cluster 
R136. Given this correspondence between our models and 
the properties of young Magellanic Cloud clusters, we are 
confident in our selection of an appropriate length scale. 

In order to allow investigation of the effects of a popu- 
lation of stellar-mass BHs on cluster evolution, we modified 
NBODY4 to allow control of the production of BHs in su- 
pernova explosions. For the present modelling this is imple- 
mented in a relatively simplistic manner; however in princi- 
ple the relevant code could be altered to cover more complex 
formation scenarios. We define three variable parameters - 
the minimum initial mass of BH progenitor stars, the masses 
of the BHs themselves, and the sizes of the natal velocity 
kicks which they receive. In each run, all stars initially more 
massive than 20 M0 produce BHs, with masses uniformly 
distributed in the range 8 < AIbh < 12 Mq so that the 
mean BH mass is 10 Mq . This range is consistent with dy- 
nami cal masses obt ained from observations of X-ray binaries 
(e.g.. lCasar"3l2006l ). We generate model clusters using the 
same random seed, so that they initially contain identical 
stellar populations. Our adopted IMF and total number of 
particles result in the formation of 198 BHs in all clusters. 

The retention fraction of BHs in a given cluster, /bh, is 
strongly dependent on the natal kicks given to the BHs at 
formation. If a kick is too strong (i.e., Ukick larger, roughly, 
than the escape velocity of the cluster, Wcsc), a BH will 
quickly cross the limiting radius of the cluster and be re- 
moved from the simulation. Under our modifications to 
NBODY4, BH kicks can vary from zero (/bh = 1) to very 
large (/bh = 0) and can be set as a constant, or selected 
randomly from a uniform distribution with specified limits, 
or a Maxwellian distribution with a specified mean. By vary- 
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ing these aspects of natal BH kicks, it is straightforward to 
control the BH retention fraction in any given model. 

Although NBODY4 allows the inclusion of primordial bi- 
nary stars in cluster models, in the present paper we investi- 
gate only models with no primordial binaries. The inclusion 
of such objects would introduce a very large new area of 
parameter space for investigation, beyond the scope of the 
time available for our simulations. Even so, any complete 
modelling of Magellanic-type clusters should undoubtably 
incorporate binary star popula tions as these are observed - 
for example. lElson et al.l (|l998l ) observed the binary fraction 
in the young massive LMC cluster NGC 1818 to be 35 ± 5 
per cent in the cluster core, decreasing to 20 ± 5 per cent in 
the outer regions. We anticipate that future simulations by 
us will investigate the effects of a binary star population on 
the results presented in this paper. 

Finally, the youngest LMC and SMC clusters typically 
have metallicities not far from the solar value - fo r exam - 
ple, the literature compilation in iMackev fc Gilmord (|2003al ) 
suggests a range -0.4 < [Fe/H] < 0.0 in the LMC. There- 
fore, for consistency, in all simulations we set our clusters 
to have solar metallicity, Z = 0.02. However, we note that 
there is a strong age-metalli city relationship present in b oth 
Magellanic Clouds (see e.g.. |Pagel fc Tautvaisienelll998h . in 
that older clusters are typically much more metal-poor than 
younger cl usters. This may h ave important implications for 
our results. iHurlev et al.l (|20Q4i l have demonstrated that dif- 
ferences in metallicity can result in some weak variation in 
the global structural and dynamical evolution of open clus- 
ters, mainly due to differences in stellar winds and mass 
loss. Furthermore, variations in metal abundance may have 
a strong effect on the numb er and mass of BH s produced in 
supernova explosions (e.g.. IZhang et al]|2007l ). We discuss 
these aspects further in Section [5] 



2.2 Primordial mass segregation 

Almost all young massive star clusters which have been ob- 
served with sufficient resolution are seen to exhibit some 
degree of mass segregation. This is true for clusters in the 
LMC (e.g., NGC 1805, NGC 1818, R136) and SMC (e.g., 
NGC 330), as well as in the Gala xy (e.g., Orion Nebu la 
Cluster, Arches, Quintuplet) (e.g.. Ide Griis et al.l l2002allbl: 
Shianni ct al. 2002 ; Malum uth fc Head 1 19941 : iBrandl et all 
1996,; Hunter et al. 1996 ; Kim et al.ll2006l ). Detailed simu- 



lation s of star cluster formation (see e.g., iBonnell fc Bate! 
I2OO6I . and references therein) are consistent with these ob- 
servations, suggesting that mass segregation in young mas- 
sive clusters may well be a product of the formation pro- 
cess, in that more massive stars are preferentially formed at 
the bottom of local potential wells where the gas density is 
greatest. 

Irrespective of whether the observed properties of young 
massive clusters are truly "primordial" , we would like to in- 
clude the possibility of very early mass segregation in our 
models in order to investigate its effects on their subsequent 
evolution. To produce initially mass segregated clusters in 
a "self-consistent" fashion (i.e., close to virial equilibrium, 
with all members having appropriate velocities) we devel- 
oped the following procedure. For a given model, we first 
generate a cluster as described in the previous Section. This 
object represents the case where there is no primordial mass 



segregation. We then implement a mass-function truncation, 
setting all stars in the cluster with masses greater than 8 Mq 
to have mass 8M0. Next, the cluster is evolved dynami- 
cally using NBODY4, but with the stellar evolution routines 
turned off. Hence the cluster begins to dynamically relax 
and mass segregate. The degree of primordial mass segrega- 
tion is then easily controlled using a single parameter - the 
length of time, Tms, for which the cluster is "pre-evolved" . 
We selected the truncation limit of 8 M0 empirically so that 
the pre-evolution can extend for a reasonable duration (a few 
hundred Myr) without the most massive stars sinking to the 
cluster centre, interacting, and ejecting each other. Once the 
desired pre-evolution time is reached, we halt the simulation, 
replace the mass-truncated stars with their original masses, 
and take the positions and velocities in the pre-evolved clus- 
ter as the initial conditions for the full run including stellar 
evolution. It is straightforward to read in the pre-evolved 
cluster using nbody4, without applying any re-scaling. 

Because we replace the mass-truncated stars with their 
original masses after the pre-evolution is complete, these 
stars (which number a few hundred in any given model) 
have slightly incorrect velocities at the beginning of the sim- 
ulation proper. However, since they almost all reside in the 
densest part of the cluster, once the full simulation begins 
these velocities change rapidly and, within the first few lo- 
cal dynamical times, become consistent with the mass dis- 
tribution in the cluster. Hence this small inconsistency has 
a negligible effect on the long-term evolution. We also note 
that during the pre-evolution a small fraction of stars escape 
from the cluster. This is usually in the form of low-mass 
stars drifting slowly across the limiting radius, after which 
they are removed from the simulation. This process is very 
gradual however, and even the clusters with the longest pre- 
evolution times (Tms = 450 Myr) always retain more than 
96 per cent of the mass of the initial non-segregated object. 
Occasionally, despite the mass-truncation of stars, a massive 
object will interact strongly with another massive object 
during the pre-evolution, and be ejected from the cluster. 
Since we are very interested in how the most massive stars 
in the cluster affect its evolution, and would like to maintain 
a high level of consistency between the BH populations of 
different model clusters, we always replace these objects at 
the end of the pre-evolution period using their positions and 
velocities from a few output times before the ejection. Since 
this is necessary for at most a handful of stars per cluster, 
the introduced inconsistencies are again negligibly small. 

The initial central densities and core radii of our pri- 
mordially mass segregated model clusters depend on the 
duration of the pre-evolution. We selected our longest pre- 
evolution times (Tms = 450 Myr) so that the resulting clus- 
ters possess properties very similar to those observed for 
R136, which is the most compact Magellanic Cloud cluster. 
These models have — 0.25 pc and log po = 4.58 (cf. Fig. 
[31). In addition to these global properties, we examined in 
detail the radial variation in mass function slope for such 
models and compared the results with those observed for 
several young Magellanic Cloud clusters. This process is de- 
scribed in detail in Section [4.21 here, we simply note the ex- 
cellent agreement between the models and the real clusters, 
as verification of the validity of our pre-evolution algorithm. 
We also ran simulations using clusters with more intermedi- 
ate pre-evolution durations (Tms = 115 and 225 Myr) - as 



8 A. D. Mackey et al. 



might be expected, these objects possess intermediate core 
radii and central densities: Tc — 0.83 pc and logpo ~ 2.70, 
and Tc. = 0.37 pc and logpo = 3.61, respectively. 



3 "OBSERVING" THE SIMULATIONS 

Since the radius-age trend is defined observationally (i.e., by 
Fig. m, a vital ingredient in our analysis is to derive mea- 
surements from the simulations which are fully consistent 
with these observations. This requirement highlights a key 
advantage in running direct, realistic A''-body models. Be- 
cause the positions, velocities, masses and luminosities of 
all stars are explicitly followed, and because we do not have 
to worry about scaling our results with A'^, we are able to 
perform simulated observations of a model cluster at each 
output time which lead to measured quantities that are di- 
rectly comparable to those obtained for the real Magellanic 
Cloud clusters. More specifically, we calculate the observa- 
tional core radius of each model cluster rather than using 
the traditional A^-body definition (see below), and further, 
we incorporate many of the subtleties of the actual HST 
measurements which have defined Fig. [1] 

Consider Fig. 31 where we have plotted the detection 
limits in the HST WFPC2 and ACS imaging from which 
Fig. [T] was constructed, against cluster age. The brighter 
limits represent saturation on the images (very bright stars, 
while recorded on the images, are generally not measured 
by photometry software), while the lower limits represent 
the approximate 50 per cent detection completeness levels 
(faint stars are not always detected above background noise 
by photometry software) . We have split the clusters into four 
age bins according to approximately constant detection lim- 
its - these are delineated on the plot with solid vertical lines. 
Within each bin, we mark the mean bright and faint detec- 
tion limits with dashed lines, and the approximate maximum 
scatter about these means with dotted lines. 

A number of things are evident from Fig. |4l First, for 
any given cluster, the observations sample only a portion 
of the range of stellar masses present in the cluster. Hence, 
the surface brightness profile, from which the structural pa- 
rameters for that cluster are measured, is based only on the 
spatial distribution of stars within this range. Second, the 
sampled range varies systematically with cluster age. This 
is due to the fact that observations of star clusters in the 
LMC and SMC are commonly aimed at targeting stars near 
the main-sequence turn-off. Consequently, the required ex- 
posure time increases with cluster age, meaning that both 
the brighter and fainter detection limits become deeper with 
age. Looking at the two oldest bins, one can also see the in- 
creased capabilities of the ACS instrument compared with 
WFPC2. While the saturation limits are comparable for all 
clusters in these two bins, the ACS-imaged objects have faint 
detection limits ~ 2 mag fainter than those of the WFPC2- 
imaged objects - indicative of the increased sensitivity and 
dynamic range of ACS over WFPC2. 

Given the above, when deriving structural measure- 
ments from our A'^-body simulations we must account for 
the fact that in the real observations for any given cluster, 
only a portion of the stellar mass function was sampled, and 
the fact that this range changes with the age of the target 
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Figure 4. Bright and faint stellar detection limits on the 
HST/WFPC2 and ACS images of LMC and SMC clusters used 
for the measurements presented in Fig. [T] Circular symbols mark 
LMC clusters, while SMC objects are trian gles. Filled symbols 
represe nt the WFPC2 imaging described in iMackev &: Gilmord 
while open symbols are the ACS imaging from Mackey 
et al. (2008, in prep.). Clusters are split into four age bins, de- 
lineated with solid vertical lines. Within each age bin, the mean 
bright and faint detection limits are marked by dashed lines, while 
the approximate maximum scatter about each mean is marked by 
a pair of dotted lines. For ease of reference, absolute magnitudes 
Mv = (—5, 0, 5) correspond to zero-age main sequence stellar 
masses of Af, ~ (45.0, 4.20, 1.06) Mq at solar metallicity. 



cluster. Since we know the details of the sampling from Fig. 
m this is simple to achieve. 

To "observe" a model cluster, we pass the A^-body data 
through a measurement pipeline essentially identical to that 
we used to obtain structural quantities for the LMC and 
SMC sample s (full details of the observ ational pipeline may 
be found in iMackev fc Gilmord [2003al ) . For the data pro- 
duced at a given output time, we first convert the lumi- 
nosity and effective temperature of each star in the clus- 
ter to UBVR I standa r d mag nitudes, using the bolometric 
corrections of Kuru c z (Tggz) supplemented with those of 
Be rgeron. Wesemael & BcauchamD (1995) for white dwarfs. 
We also convert the position and velocity of each star to 
physical units using the appropriate length and velocity scale 
factors (see Section[2TT} ■ With this complete, we next impose 
the bright and faint detection limits appropriate to the out- 
put time (these are the dashed mean limits in Fig. [4]). This 
leaves an ensemble of stars with which t o construct a surface 
brightn ess profile, which we do following IMackev fc Gilmord 
(|2003al lbl). We project the three-dimensional position of each 
star onto a plane, construct annuli of a given width about the 
cluster centre, and calculate the surface brightness in each 
annulus. For consistency with the observational pipeline, we 
use a variety of annulus widths so that both the bright inner 
core and the fainter outer regions of the cluster are well mea- 
sured. Next, we account for the fact that both the WFPC2 



Core expansion in massive star clusters 9 



and ACS cameras have fields-of-view which are considerably 
smaller than the area on the sky filled by a Magellanic Cloud 
cluster, which might typically have rt ~ 40 — 50 pc (i.e. , 
rt ~ 160-200" a t the LMC distance) (see e.g.. lMateolll987l : 
lOlsen et al.lll998l ). This results in surface brightness profiles 

G enerally being truncated b eyond projected radii rp ~ 25 pc 
Mackev fc Gilmore|[2003al lbfi. After imposing this limit, we 
finally fit an EFF model of the form of Eq. [T] to the result- 
ing surface brightness profile, and from this model derive 
the structural parameters - in particular the core radius, r^, 
and the power-law slope at large radii, 7. To reduce noise we 
repeat this process for three orthogonal planar projections 
at each output time and average the results. 

It is worth noting the difference between the quantity 
and the 'core radius' usually defined in A'^-b ody simulations. 
This has been discussed in some detail by IWilkinson et al.l 
however, in the interests of clarity we re-iterate a 
few of the most salient points. Traditionally, observers, the- 
orists, and numericists have employed difi'erent interpreta- 
tions of the 'core radius'. That for observers is as defined 
above (Eq. [2| , as the projected radius at which the surface 
brightness (or density) has dropped to half the central value. 
Theoretically defined, the core radius is the natural scale- 
length of the model under consideration - for example, in 
EFF models a is the scale-length. Eq. [2] provides a general 
relation between a and Vc. It should be noted however, that 
as a cluster evolves, the EFF parameters are not static, and 
therefore the ratio between a and Tc is variable with time. 

In A'-body simulations the numerically calculated 'core 
radius' is more correctly termed the density radius, rd- The 
im plementation in NBODY4 is based on a quantity described 
by ICasertano fc Hud l|l985D . so that Td is defined as the 
density-weighted average of the distance of eac h star from 
the density centre of the cluster l|Aarsethll200ll ). The local 
density at each star is computed from the mass within the 
sphere containing the six nearest neighbours. As noted by 
IWilkinson et all l|2003l ), there is no general relationship be- 
tween Vd and rc, and in fact the behaviour of and ra may 
be quite different throughout a simulation. 

As a final remark, we briefly consider the appropriate- 
ness of fltting a power-law proflle (Eq. [T|, which formally 
has no outer limit, to a simulated cluster evolving in a tidal 
fleld. There are two reasons why this is acceptable. First, 
because of the radial truncation imposed to mimic the field- 
of-view limitations of the WFPC2 and ACS cameras, our 
derived surface brightness profil es do not reach as far as the 
cluster tidal radius. Following M ackev fc Gilmor3 (|2003al lbh. 
it is therefore legitimate to fit EFF models to these observed 
profiles, even when a cluster is dynamically old enough to 
exhibit a tidal truncation - in the interests of obtaining mea- 
surements of 7 which are, like those for Vc, directly compa- 
rable to the real observations, we choose to employ the same 
methodology. Even without the truncation of our radial pro- 
files, an EFF model would still have been the most appro- 
priate choice. This is due to the treatment of stellar escapers 
in NBODY4, as discussed in Section [2.11 While the tidal ra- 
dius rt is estimated from Eq. |4l stars are not removed from 
the simulation until they reach 2rt. Hence they are free to 
populate the region rt < r < 2rt, and there is no truncation 
in the density profile at (or near) rt, even for dynamically 
old clusters. 



4 SIMULATIONS AND RESULTS 

The properties of our A-body runs are listed in Table [T] 
Our main set of models are labelled Runs 1-4. These cover 
the extremes of the parameter space we are interested in 
investigating, spanned by BH retention fractions /bh = 
and /bh ~ 1, and the pre-evolution durations Tms ~ 
Myr (i.e., no primordial mass segregation) and Tms ~ 450 
Myr (strong primordial mass segregation, matching that ob- 
served in young LMC and SMC objects). These runs are 
therefore expected to represent the extremes of cluster evo- 
lution induced by variation of the BH retention fraction and 
the degree of primordial mass segregation. The global prop- 
erties of thes e four Runs have al ready been presented in a 
short Letter (|Mackev et al.l [20071 ): in the present paper we 
examine their evolution in considerably more detail. 

In addition to our four primary runs, we performed sev- 
eral additional simulations in order to sample the parameter 
space more completely, and in particular verify that models 
with intermediate values of /bh and Tms exhibited evolution 
intermediate between that displayed by Runs 1-4. To this 
end. Runs 4a and 4b explore the effects of primordial mass 
segregation in more detail, while Run 5 highlights the effects 
of natal kicks on BH retention and the subsequent cluster 
evolution. Finally, Run 6 is used to address the question of 
whether we can reproduce the cluster evolution induced by a 
significant BH population by retaining neutron stars (NSs) 
instead of the BHs. 

For each run, we measured the initial cluster mass, cen- 
tral density, and observed structural parameters rc and 7 - 
these are all listed in Table[T] It is important to re-emphasize 
how closely these correspond to the observed quantities for 
the youngest massive clusters in the Magellanic Clouds. This 
can be seen explicitly by comparing the values listed in Ta- 
ble [T] with the plots in Fig. O The model clusters with no 
primordial mass segregation have rc ~ 1.9 pc, 7 ~ 3.0, and 
log po ~ 2.3. These clusters therefore appear very similar 
to a number of Magellanic Cloud clusters with ages of ~ 20 
Myr. In contrast, the heavily mass segregated model clusters 
have much smaller cores and higher central densities, with 
rc ~ 0.25 pc and logpo ~ 4.6. They also have flatter power- 
law fall-offs, with 7 ~ 2.3. In this respect, they strongly re- 
semble the very compact massive young LMC cluster R136, 
which has an age of ~ 3 Myr. The total masses of all mod- 
els are very similar, in the range 4.728 log Mtot ^ 4.746. 
The variation is due to the pre-evolution procedure used to 
develop the mass segregated initial conditions, as described 
in Section [2.21 From comparison with Fig. (3] it is clear that 
our A-body clusters have masses typical of the youngest 
clusters in the observed sample. We also note that the "ob- 
served" integrated colours of our models at early times are 
consistent with measurements for young Magellanic Cloud 
clusters from the literature - for example, th e integrated 
{B — V) colours compiled bv lBica et al.l l| 19961 '). 

Given the close correspondence between the properties 
of our model clusters and those observed for young LMC and 
SMC objects, we are confident that our A-body simulations 
are directly modelling the evolution of massive Magellanic 
Cloud clusters. 

Output data was produced for each Run at intervals 
of At — 1.5 Myr for r <= 100 Myr, and at intervals of 
At = 15 Myr for r > 100 Myr. It is worth noting that no 
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Table 1. Details of A'^-body runs and initial conditions. Each cluster begins with A'^o stars with masses summing to Mtot, and 
initial central density po- Initial cluster structure is "observed" to obtain and 7. Each model is evolved until Tmax- 



Name 


Nq 


log Mtot 


log po 




7 


Initial MSeg 


BH Retention 








(Mq) 


(Mq pc-3) 


(pc) 




(2ms) 


(/bh) 


(Myr)l 


Run 1 


100 881 


4.746 


2.31 


1.90 ±0.09 


2.96 ±0.17 


None 


0.0 


19 987 


Run 2 


100 881 


4.746 


2.31 


1.90 ±0.09 


2.96 ±0.17 


None 


1.0 


10 668 


Run 3 


95 315 


4.728 


4.58 


0.25 ±0.04 


2.33 ±0.10 


450 Myr 


0.0 


11274 


Run 4 


95 315 


4.728 


4.58 


0.25 ±0.04 


2.33 ±0.10 


450 Myr 


1.0 


10 000 


Run 4a 


98 605 


4.738 


2.70 


0.83 ±0.07 


2.45 ±0.14 


115 Myr 


1.0 


4274 


Run 4b 


97 209 


4.733 


3.61 


0.37 ±0.05 


2.34 ±0.10 


225 Myr 


1.0 


4 457 


Run 5 


95 315 


4.728 


4.58 


0.25 ±0.04 


2.33 ±0.10 


450 Myr 


0.5 


10 059 


Run 6 


100 881 


4.746 


2.31 


1.90 ±0.09 


2.96 ±0.17 


None 


0.0, NS2 


19 987 



^ As described in Section |4] no special significance should be attached to the listed values of rma: 
^ Run 6 is identical to Run 1, except with natal neutron star kicks set to zero so that /ns = 1-0. 



special significance should be attached to the listed values 
of Tmax in Table [T] The main criterion for our primary Runs 
(Runs 1-6, excluding Runs 4a and 4b) was that Tmax be 
larger than ~ 10 Gyr, to approximate the ages of the oldest 
Magellanic Cloud globular clusters. The listed Tmax simply 
represent the most convenient termination points beyond 
this time. For interest's sake, Runs 1 and 6 were evolved for 
significantly longer periods (rmax = 20 Gyr) than the other 
models, so that the clusters passed through the core-collapse 
phase. In contrast. Runs 4a and 4b were evolved only as 
long as necessary (i.e., just long enough for the effects of 
intermediate values of Tms to become evident), to save on 
computation time. 



4.1 Runs 1 and 2: No mass segregation 

We first consider the pair of simulations labelled Run 1 and 
Run 2. Neither of these two model clusters have primordial 
mass segregation, and both start with identical initial con- 
ditions, to the extent that they share the same random seed. 
The sole difference between them is that in Run 1 the natal 
BH kicks are set to be Vkick ~ 200 kms~^, whereas in Run 2 
they are set to be zero. Thus, every BH formed in a super- 
nova explosion in Run 1 is provided with a sufficiently large 
random velocity that it very rapidly escapes from the clus- 
ter, so the retention fraction is /bh = 0. Conversely, in Run 
2 all 198 BHs are retained in the cluster and the retention 
fraction is /bh = 1- The purpose of these runs is twofold. 
First, Run 1 lets us consider the long-term evolution of our 
simplest cluster set-up - no primordial mass segregation, and 
zero BH retention. This model therefore constitutes a con- 
trol run against which the evolution of all our other models 
may be compared. Second, by making such a comparison. 
Run 2 lets us isolate the effects of a population of stellar- 
mass black holes on the structural and dynamical evolution 
of a massive star cluster. 

The progress of Runs 1 and 2 across the radius-age plane 
is displayed in Fig. [5] (left panel) . Also shown is the evolu- 
tion of these two runs in the 7-age plane (right panel). First 
consider Run 1, which behaves exactly as expected for a 
classical massive stellar cluster. At very early times, extend- 
ing to roughly r ~ 100 Myr, there is a period of severe 



mass loss due to the rapid evolution of the most massive 
stars in the cluster. By r ~ 100 Myr, approximately 25 per 
cent of the initial cluster mass has been lost. The 198 BHs 
are formed in supernova explosions between 3.5-10 Myr and, 
since they are born with Ukick ~ 200 kms~^, all are immedi- 
ately ejected from the cluster. From Fig.O it is clear that the 
violent relaxation experienced by the cluster when r < 100 
Myr is not refiected in its core- radius evolution, presumably 
because the mass loss is distributed evenly throughout the 
cluster. Similarly, there is no evidence of the violent relax- 
ation phase in the evolution of 7. 

As the cluster grows older, the rate of mass loss de- 
creases and the cluster settles into a quasi-equilibrium state, 
where dynamical evolution is dominated by two-body re- 
laxation processes. The median relaxation time for this 
N = 10° star cluster is give n by Uh « 1.9x10^^ Af^Y^mjV^/^ 
(|Binnev fc Tremaind [19871 ') where m* is the typical stellar 
mass and is the 3-dimensional radius containing 0.5Mci. 
At r = 100 Myr, when the rapid early mass loss is mostly 
complete, Md « 43 500 M©, m* « 0.45 M© and rh ^ 8 pc, 
so that trh ~ 2 Gyr. Mass segregation develops in the cluster 
on roughly this time-scale: this is evident in Fig.[5]as a grad- 
ual contraction in Vc as the most luminous stars in the mag- 
nitude range used to measure the structural parameters (cf. 
Fig. |4]) sink towards the cluster centre. As two-body relax- 
ation proceeds and mass segregation becomes more promi- 
nent, the core radius steadily shrinks with time. The power- 
law fall-off, 7, slowly becomes steeper during this phase; 
however as the core becomes increasingly more compact, so 
7 becomes increasingly flatter after r ~ 10 Gyr. 

Eventually, after many Gyr of evolution. Run 1 en- 
ters the core-collapse phase. The point of greatest collapse 
(smallest Vc) occurs at r « 17.4 Gyr, when the central 
mass density reaches logpo « 4.5 - a value commensurate 
with those inferred for NGC 2005 and 2 019, the most hkely 



core c o llapsed clusters in the LMC (e. g., Mackey fc Gilmorg 



l2003a|; iMcLaughlin" fc van der MarellBoOSh . The point of 
greatest collapse coincides with a spate of binary star for- 
mation in the core - by r = 17.5 Gyr there are seven newly- 
formed binary stars. Subsequently, up until the end of the 
simulation at r = 20 Gyr, there is no significant change in 
the observed value of Vc- Defining the cluster age in terms 



Core expansion in massive star clusters 11 




6789 10 6789 10 

logio(^/Myr) logjo(T/Myr) 



Figure 5. Structural evolution of Af-body Runs 1 and 2. Neither has primordial mass segregation; the only difference between them 
is the BH retention fraction (/bh = and 1, respectively). Left panel: Evolution of rc, observed as described in Section [3] Run 

1 evolves exactly as expected for a classical massive star cluster, with the main trend being a slow contraction in Tc as the system 
relaxes dynamically and moves towards core collapse. In stark contrast, Run 2 evolves very similarly to a point, after which strong 
expansion in the core radius is observed. The presence of 198 stellar-mass black holes in this cluster thus leads to strikingly different 
core radius evolution. Right panel: Evolution of the power-law fall-off, 7, again observed as described in Section|3] As with rc. Run 2 
evolves identically to Run 1 until the BH population becomes dynamically active, after which the evolution strongly diverges with Run 

2 devel oping a steep fall-off in its outer regions. In this panel, only data points from the WFPC2 observations of Mackcy & Gilmor^ 
ll2003al lbl) are plotted. Measurements of 7 from the recent ACS observations of Mackey et al. (2008, in prop.) are not yet finalised. 



of an integrated median relaxation time, which is necessary 
because trh is a constantly evolving quantity, we find that 
at r = 17.4 Gyr, 8.37 trh has elapsed. 

Fig. |6] illustrates the evolution of Run 1 in more detail. 
In the top panel a series of Lagrangian radii are plotted - 
here, we define Rx% to be the 3-dimensional radius contain- 
ing X per cent of the stellar mass in the cluster - that is, 
excluding BHs. This exclusion is not important for Run 1, 
since all BHs are gone from the cluster by ~ 15 Myr; how- 
ever it is crucial for examining the evolution of the stellar 
component of Runs in which /bh > 0. Unlike rc, the inner- 
most Lagrangian radii in Run 1 do show an increase in size in 
reaction to the early mass-loss phase; however, this increase 
is only very modest. In addition, the innermost Lagrangian 
radii do not show any sign of contraction until much later 
than does r^ - this is an indication of the luminosity (and 
hence mass) weighting inherent in the calculation of rc. The 
half-mass radius of the cluster shows only a small amount of 
variation throughout its evolution. The outer radii also show 
only very gradual evolution. The main feature is an expan- 
sion in the 90% radius during the mass-loss phase. This is 
due to stars in the very outer regions of the cluster drift- 
ing beyond rt as the cluster rapidly loses mass. Eventually 
these objects are removed from the simulation (once they 
get beyond 2rt) and the 90% radius slowly contracts. This 
contraction continues as the cluster slowly loses mass for the 
rest of its lifetime, and rt gradually shrinks accordingly. 

In the lower panel of Fig. [S] the mean stellar mass in 
shells encompassed by selected Lagrangian radii is plotted. 
This plot therefore shows the development of mass segrega- 



tion in Run 1. This process is inhibited by the early violent 
relaxation phase, and there is only a very small degree of seg- 
regation present in the cluster's central regions by r = 100 
Myr. Subsequently however, the stratification becomes very 
well established. As expected, this occurs more rapidly in 
the central regions of the cluster, where the relaxation time 
is shortest. By the time the core-collapse phase is reached, 
there is a large degree of mass segregation present in the 
model cluster. One notable feature, exhibited by both the 
outermost shell and the full cluster mean, is an increase in 
the mean stellar mass after r ~ 10 Gyr. This is due to the 
preferential removal of low-mass cluster stars by the external 
tidal field. These stars typically reside in the outer cluster 
regions at late times, and are hence far more susceptible to 
tidal effects than are the more massive objects which inhabit 
the cluster core. At late times, stellar evolution has all but 
slowed to a halt so that the tidal stripping of low-mass stars 
has a significant effect on the mean stellar mass. 

Now consider Run 2, in which the BH kicks Wkick = 
Okms^^ so that /bh = 1- From inspection of Fig. (5] it is 
clear that the evolution of both rc and 7 is indistinguishable 
from that seen in Run 1 up to log r ~ 8.8, after which strong 
expansion of rc is observed for Run 2, in conjunction with 
a significant steepening in 7. A careful comparison between 
the two models reveals that this divergence begins at r ~ 
650 Myr. The expansion of rc in Run 2 continues for the 
remainder of the simulation, until Tmax = 10.67 Gyr. Since 
Runs 1 and 2 are identical apart from the kicks imparted to 
the BHs on their formation, the strongly different evolution 
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Figure 6. Evolution of various Lagrangian radii (top panel) and 
the mean stellar mass in the shells encompassed by selected La- 
grangian radii (lower panel) for Run 1. The radii displayed in 
the top panel are, from inner to outer, the 1%, 5%, 10%, 30%, 
50% = Th (dashed line), 70%, 80%, and 90% radii. In the lower 
panel the shells are defined by: r ^ Ri%^ Rl% < r ^5%i 
^5% < r ^ Rio%' RlO% < r ^ ^30%. and Rjq% < r ^ Rgo% 
(these are listed in order from the upper to lower solid lines at 
the right hand side of the panel). The dashed line is the mean 
mass for all stars in the cluster. 



observed for these two models must be due to the presence 
of the 198 BHs in Run 2. 

The properties of this BH population as a function of 
time are illustrated in Fig. [T] As in Run 1, by r ~ 100 Myr, 
the most violent phase of stellar evolution is essentially com- 
plete. At this time, the BHs (of typical mass ttibh = 10 Mq) 
are already significantly more massive than any other cluster 
members (of typical mass m, ~ 0.45 Mq), and are hence be- 
ginning to sink to the cluster centre via dynamical friction, 
on a time-scale of ~ (m»/mBH) trh ~ 90 Myr. 

This is evident from panels b and c in Fig. [T] Panel b 
shows the number of BHs within the shells encompassed by 
the stellar Lagrangian radii r ^ Ri%, R\% < r ^ RsVo, and 
r > Rio%. The evolution of these Lagrangian radii them- 
selves may be seen in Fig. |8l which is discussed in more 
detail below. Panel c shows the evolution of the black hole 
Lagrangian radii -Bio%, ^25%, -850%, and -875%, where, by 
analogy with the stellar Lagrangian radii, is the 3- 

dimensional radius containing x per cent of the BH mass 
in the cluster. 

By 200 Myr, the mass density of BHs at the centre of 
the cluster is already roughly equal to that of the stars, and 
by 400 Myr it is about three times larger. Shortly after, this 
central B H subsystem becomes unstable to further contrac- 
tion (see ISpitzerlll987l . Eq. 3-55) and decouples from the 
stellar core in a runaway gravothermal collapse, leading to a 
very rapidly increasing central BH density - by 490 Myr, the 
central density of the BH subsystem is ~ 80 times that of 




T (Myr) 

Figure 7. Properties of the BH population in Run 2 as a function 
of time: (a) the number of single BHs (upper line) and binary BHs 
(lower line) in the cluster; (b) the number of BHs within the shells 
encompassed by the steHar Lagrangian radii (cf. Fig. |6ll r ^ Ri%, 
R\% < "T ^5% I and r > Riq% (the upper, middle and lower 
lines, respectively, at the right of the plot); (c) the black hole 10%, 
25%, 50% and 75% Lagrangian radii (respectively, the innermost 
to outermost lines); (d) the cumulative numbers of escaped single 
BHs (upper) and binary BHs (lower), along with fits of the form 
Ne = Ao + A\T — yl2TlogT (dashed lines); and (e) the radial 
positions of three typical BHs. The vertical dotted line indicates 
T = 650 Myr, the approximate time when core expansion begins. 
The evolution of Tc is plotted (dot-dashed line) in panel (e). Note 
the different axis scales on either side of panels (a) and (d). 



the stars. This is sufficiently dense for the creation of stable 
BH binaries in three-body interactions to be initiated - the 
first such object is formed at ~ 510 Myr, and by ~ 650 Myr 
there are several (see Fig. [7^). At this point, the collapse of 
the BH subsystem is halted: the BH Lagrangian radii cease 
their inward movement and become roughly constant, while 
the number of BHs within the inner stellar Lagrangian radii 
also level off. It is at this time that the evolution of the ob- 
servational structural parameters Tc and 7 in Run 2 begins 
to strongly deviate from that in Run 1. 

As noted above, prior to this point the evolution of Run 
2 is observationally identical to that of Run 1. Neither the 
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retention of the BH population at r « 10 Myr, nor the 
subsequent orbital decay of these objects and the result- 
ing formation of a compact central BH subsystem leads to 
differential evolutio n of rc- This appears at odds with the 
models presented bv lMerritt et alj (|2004l '). who investigated 
the possibility that the radius-age trend results from the for- 
mation of cores in primordially cusped star clusters due to 
the sinking and central accumulation of massive stellar rem- 
nants. We attribute our differing results to the much higher 
degree of cent r al ma ss concentration in the cusped models of 
iMerritt et alJ (|2004h . which thereby respond more strongly 
and more rapidly to the perturbations induced by sinking 
remnan ts than does o u r initi ally cored, relatively low density 
Run 2. iMerritt et all (|2004 ) also mentioned the possibility 
of additional cluster expansion due to the subsequent evolu- 
tion of the BH population, once the central subsystem had 
formed. As discussed below, all of the expansion observed in 
our models is the result of such processes. 

The number of stable BH binaries in Run 2 peaks at 5 
at r ~ 890 Myr. After this point, there are — 5 BH bina- 
ries at any given time (Fig. [TK). Once formed, a BH binary 
undergoes superelastic collisions with other, usually single, 
BHs in the central core (although BH binaries do also occa- 
sionally collide with each other). On average, as BH binaries 
participate in such interactions they become "harder" (more 
tightly gravitationally bound) , with the released bin ding en- 
ergy being carried off by the interacting BHs (e.g.. iHeggid 
ll975l : lHeggie fc Hutll2003l '). In each such interaction, the bi- 
nary BH also has a recoil velocity imparted to it, the magni- 
tude of which is dependent on how energetic the interaction 
has been. Together, these processes result in the scattering of 
BHs outside r^, often into the cluster halo. As a given binary 
becomes increasingly tightly bound, so too can the collisions 
in which it is involved become increasingly energetic, such 
that an interacting BH carries off sufficient kinetic energy 
that it escapes from the cluster altogether. Eventually the 
BH binary is sufficiently hard that the recoil velocity it re- 
ceives during a collision is larger than the cluster escape 
velocity, and the binary escapes as well. Hence, interactions 
in the central compact BH subsystem also result in the ejec- 
tion of BHs from the cluster. For clarity we will retain the 
italicised terminology {scattering and ejection) henceforth. 

These processes are evident in Fig. [7^, which shows the 
movement of three typical BHs during Run 2. Each of the 
three is born well outside Tc, but all sink to the central core 
via dynamical friction, as described above. Two are already 
present there by the time the first BH binaries are formed. 
All three of the BHs are frequently scattered to Tc (dot- 
dashed line) during their evolution in the cluster, and at 
least once each into the cluster halo. One is ejected from 
the cluster ai t — 7900 Myr due to a strong interaction 
in the core. Another becomes a member of a BH binary at 
r = 6200 Myr, and subsequently undergoes four strong in- 
teractions (including one in which its partner is exchanged) , 
with increasing recoil velocity each time until this is suffi- 
cient for ejection at r = 8200 Myr. 

Considering Fig. [7};, it is clear that at any given time 
there are always a handful of BHs outside the 10 per cent 
stellar Lagrangian radius. This is an indication of the on- 
going scattering of BHs to outside ~ Vc, since any ejected 
BHs tend to escape the cluster quite rapidly. As is evident 
from Fig. [7^, a scattered BH gradually sinks back into the 
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Figure 8. Evolution of various Lagrangian radii (top panel) and 
the mean stellar mass in the shells encompassed by selected La- 
grangian radii (lower panel) for Run 2. The radii displayed in 
the top panel are, from inner to outer, the 1%, 5%, 10%, 30%, 
50% = (dashed line), 70%, 80%, and 90% radii. In the lower 
panel the shells are defined by: r ^ Ri%, Ri% < r ^5%i 
^5% < r ^ Rio%^ RiO% < r ^ ^30%, and R^o% < r ^ _R8o% 
(these are listed in order from the upper to lower solid lines at 
the right hand side of the panel). The dashed line is the mean 
mass for all stars in the cluster. The evolution of Run 1 is marked 
by dotted lines, for comparison (note that the abscissa does not 
extend to such late times as are plotted in Fig.|6ll. 



cluster centre via dynamical friction, thus transferring its 
newly-gained energy to the stellar component of the clus- 
ter. Most is deposited within Vc, where the stellar density 
is greatest. The ejection of BHs also transfers energy to the 
cluster, since a mass m escaping from a cluster potential well 
of depth |$| does work m|$| on the cluster. This mechanism 
is particularly effective in heating the stellar core, since BHs 
are ejected from the very centre of the cluster, and the en- 
ergy contributed to each part of the cluster is proportional 
to the contribut ion which that part makes to the central po- 
tential (see e.g.. lHeggie fc Hutll2003l ). In addition, here mBH 
is significantly larger than m«. 

Together these two processes heat the stellar core of the 
cluster, resulting in significant core expansion. This becomes 
evident observationally at t ~ 650 Myr, and continues as 
long as the BH population is dynamically active - in the case 
of Run 2, the simulation was halted before the expansion 
ceased. From Fig. O the size of rc is roughly proportional 
to log r, consistent with the shape of the upper envelope of 
the observed cluster distribution. However, in this A''-body 
model the expansion begins too late for the evolution to 
trace the upper envelope exactly; rather, it runs parallel. 

The evolution of the stellar Lagrangian radii in Run 2 
is illustrated in Fig. |S1 along with the evolution of the mean 
stellar mass in the same Lagrangian shells examined earlier 
for Run 1. The progress of Run 1 is also marked on Fig. [8] 
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for comparative purposes (dotted lines) . As noted above, the 
initial infall and accumulation of BHs in the cluster centre 
does not cause any differential expansion of Run 2 over Run 
1 at any radii. It is only after BH binaries are formed and 
the BH population becomes dynamically active that expan- 
sion occurs in Run 2. This expansion is evident at all radii, 
although the magnitude is greatest in the central regions 
of the cluster. None of the Lagrangian radii expand by as 
great a factor over the simulation as does rc. The explana- 
tion for this can be seen in the lower panel of Fig. |8] - the 
development of mass segregation amongst the stellar com- 
ponent in Run 2 is severely inhibited by the activity of the 
BIf population, compared to Run 1. This results in a larger 
apparent expansion in Vc than in the innermost stellar La- 
grangian radii because of the luminosity weighting inherent 
in the measurement of rc. 

The process of mass segregation in Run 2 is only sup- 
pressed after the BH population becomes dynamically ac- 
tive. Up until this point, segregation has been proceeding 
just as in Run 1; however, after r « 650 Myr, no further 
stratification occurs. Stellar evolution subsequently reduces 
the mean mass in each Lagrangian shell with time. This 
interpretation is consistent with the cluster expansion pro- 
cesses due to BH scattering and ejection which were de- 
scribed above. In particular, the repeated BH scattering- 
sinking cycles constantly stir up the stellar component of 
the cluster and hence hinder the development of mass segre- 
gation, particularly in the inner cluster regions. The strat- 
ification which occurs before the BH population becomes 
dynamically active is not reversed however - there is still 
clearly a mean-mass gradient from the inner to the outer 
regions of the cluster at all times. 

A useful quantity for examining the evolution of 
the cluster structures in Runs 1 and 2 is the ratio of 
the core radius to half-light (or half-mass) radius (e.g.. 



Vcspcrini fc Chcrnoflj [l993 : iTrenti et al.]|2007l : iHeggie et all 
2007 ; Hurlov..2007i ). In regards to the latter of these radii, the 



relevant observational parameter is the projected radius con- 
taining half the cluster light {r^.i)- This is straightforward to 
calculate for the EPF family of models. The enclosed lumi- 
nosity as a function of projected radius rp may be obtained 
by integrating Eg. I A2I within a cylinder of radius rp alon g 
the line of sight (e.g., Eq. 11 in lMackev fc Gilmorell2003ah : 



L{rp) 



27r^o 

7-2 



(r-p + a^) 



-(7-2)/2 



(5) 



When rp = rh,i the enclosed luminosity is half of the total 
luminosity, i.e. L{rh,i)/ Ltot ~ 1/2. Substituting this into Eq. 
[5] and rearranging the result leads to an expression for r^.i'- 



log(rL+a') = ^logf^ 



Ltot(2 - 7) ^ ^2 

■ilVflo 



(6) 



Projected half-light radii, along with the ratios rc/rh,i may 
h ence be calculated for the L MC and SMC cluster samples 
of lMackev fc Gilmord (|2003al lbh using their best-fitting EFF 
models and total luminosity estimates. Directly comparable 
quantities may also be calculated at each output time for 
our A^-body Runs using the "observed" EFF models. In this 
procedure, for the purposes of direct comparison we do not 
use Ltot as calculated by the A^-body code, but rather the 
total luminosity enclo sed within some limi t ing ob servational 
radius, as specified in iMackev fc Gilmord l|2003al ). 
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Figure 9. Evolution of the ratio of core radius to projected half- 
light radius rc/r^j for Af-body Runs 1 and 2, compared with 
measurements for LMC and SMC clusters. Note that the mea- 
sured ratios for the oldest, most compact LMC clusters are upper 
limits, reflecting the upper limits to the core radius measurements 
for these clusters (cf. Fig. [TJ. 



The evolution of rc/r^.i for Runs 1 and 2, compared 
with the measurements for LMC and SMC clusters, may 
be seen in Fig. [9] For much of Run 1, this ratio is a sta- 
ble quantity at rc/rh,i ~ 0.45. As this model enters core 
collapse, however, the ratio shrinks to become very small. 
This is very similar beh aviour t o that observed by previous 
authors - in particular iHurlevI |20o3), who measured the 
evolution of an identical (observationally defined) quantity 
in his large A^-body models. Very different behaviour is ob- 
served for Run 2, however. As soon as the BH population 
in this model becomes active and core expansion begins, 
rc/rh,i begins to steadily increase. This presumably refiects 
the increased heating efficiency of the BH population within 
the stellar core, as compared with the heating efficiency at 
larger radii in the cluster (cf. Fig. [8}. By the end of Run 
2 rc/rh.i ~ 0.8, matching the values observed for several 
of the most extended Magellanic Cloud clus ters. These ob - 
servations are consistent with the results of *H urle"vl |20o3), 
who found that even the presence of one BH-BH binary can 
prevent the expected decrease in rc/rh,i - in our models the 
presence of many BHs results in a significant increase in this 
ratio. It has been suggested that a cluster with a large value 
of rc/rh,i may harbour a central IM BH (see the extensive 
discussion presented by lHurlevl[2007l ): however, our Run 2 
clearly demonstrates that the presence of a population of 
stellar-mass BHs can also lead to large values of this ratio. 

Returning to Fig. [H] one significant point of note is that 
although the spatial distribution of stellar mass is quite dif- 
ferent in Run 2 compared to Run 1, the overall mean stellar 
mass in Run 2 (dashed line) remains almost exactly the same 
as that in Run 1 throughout the simulation. Given that the 
initial stellar populations in the two models were identical, 
this indicates that the typical mass of a star escaping across 
the tidal radius in the two runs is very similar. Calculat- 
ing the mean mass of all escaping stars in Run 1 between 
r — 100 Myr (when the early violent relaxation is essentially 
complete) and r — 10667 Myr (when Run 2 is terminated) 
reveals a value of 0.328 M0, while the same calculation for 
Run 2 results in 0.332 Mq. These two values are indistin- 
guishable, which is remarkable given the strong divergence 
in the structural evolution of the two clusters. Inspection 
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of the distribution of velocities with which stars escape be- 
tween r = 100 — 10667 Myr in each simulation reveals these 
also to be indistinguishable. Together these results imply 
that both models lose stars solely due to relaxation pro- 
cesses. There is only a tiny group of ~ 20 relatively high 
velocity stellar escapers in Run 2 (i.e., which have an escape 
velocity greater than that of the fastest escaper in Run 1) 
out of a total of more than 55 000 stellar escapers, indicat- 
ing that stars interact closely with BH binaries only very 
rarely. Heating of the stellar component via close interac- 
tions between stars and BH binaries is therefore negligible 
- the hardening of BH binaries is driven solely through in- 
teractions with other BHs in the central subsystem. 

It is also enlightening to consider the properties of the 
escaping BHs in Run 2. The cumulative number of escaped 
single and binary BHs is plotted in Fig. [7ji. The approxi- 
mate time at which core expansion begins, r = 650 Myr, 
is marked with a vertical dotted line. Some single BHs es- 
cape before this point - these are BHs which are formed in 
the outer regions of the cluster and drift across the tidal 
boundary due to the early violent fluctuations in the clus- 
ter's gravitational potential. After r = 650 Myr, once BHs 
begin to be ejected solely due to interactions in the central 
subsystem, it is clear that the cumulative numbers A^e of es- 
caping single and binary BHs increase more slowly at later 
times - that is, that the escape rates decrease with time. 
Hypothesising that the time derivatives of these rates vary 
as — 1/t (i.e., that dNc/dr oc — logr) suggests a fit of the 
form Nei^r) = Ao + A\t — A2T log r to the cumulative distri- 
butions, where the Ai are coefficients derived in the fitting 
process. Best-fit curves of this form are also plotted in Fig. 
[TJi (dashed lines) . Clearly these are excellent matches to the 
observed cumulative distributions, indicating that the rates 
of single and binary BH escape do indeed both have time 
derivatives which vary as — 1/t. 

The BH escape rates decrease with time because the 
density of the central BH subsystem is also decreasing with 
time - this is evident from Fig. [T];, which shows that the in- 
ner BH Lagrangian radii follow a generally increasing trend 
throughout the majority of the simulation. The typical num- 
ber of BHs in the central subsystem falls with time be- 
cause of BH ejections (Fig. [T^j), and these ejections also 
heat the BH core. Simultaneously, the stellar component of 
the cluster is becoming more extended, meaning that the 
gravitational potential at the centre of the cluster due to 
this component is becoming increasingly shallow. Together 
these processes lead to the density of the central BH sub- 
system decreasing, on average, with time. The mean BH-BH 
encounter rate also therefore decreases with time, meaning 
that the BH binary hardening rate decreases, as does the 
BH ejection rate, as observed. 

The decreasing BH binary hardening rate also means 
that the BH scattering rate decreases with time. Together 
with the slowing BH ejection rate, this means that the stel- 
lar core is also less efficiently heated with time. This is re- 
flected in the roughly logarithmic dependence of rc on r. Be- 
cause the BH scattering and ejection rates decrease through- 
out the lifetime of Run 2, by the end of the simulation at 
Tmax = 10.67 Gyr, there is still a sizeable population of 
65 single BHs and 2 binary BHs remaining in the cluster. 
This contrasts strongly with the results from early, more 
analytic, studies of the evolution of BH subsystems in glob- 
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Figure 10. Separations and eccentricities of the ejected BH bi- 
naries in Run 2 as a function of cluster age. Eccentricity is repre- 
sented by point style: BH binaries with e ^ 0.8 are asterisks, those 
with 0.8 < e ^ 0.95 are open circles, while those with e > 0.95 
are filled circles. The asterisk marked with 'Tr' is the innermost 
binary in the one ejected triple BH system (see text). 



ular clusters, which predicted depletion of any BH popu- 
lations on timescales much less than the cluster lifetimes 
llKulkarni. Hut fc McMiUanI 1 19931 : ISigurdsson fc Hernguistl 
Il993tl . The fact that the BH encounter rate decreases due 
to the interplay between the stellar component of the clus- 
ter and the BH population, as seen in our detailed numerical 
modelling, prolongs the life of the BH subsystem in a massive 
star cluster for much longer than previously appreciated. 

The properties of the ejected BH binaries in Run 2 com- 
plete the picture of BH evolution in this model. Over the 
course of the simulation, 15 BH binaries are ejected. Their 
separations [ab) and eccentricities (e) are displayed in Fig. 
1101 as a function of cluster age. The hardest binaries are 
clearly ejected at the earliest times. This is when the cluster 
escape velocity (wcsc) is largest - that is, when the binaries 
can be hardened to the greatest extent before the recoil ve- 
locity imparted during close interactions ejects them from 
the cluster. Typical separations are a;, « 4 — 6 AU for bi- 
naries ejected when r < 2.5 Gyr. As the cluster expands 
and loses mass, Vasc decreases and BH binaries are ejected 
before becoming this hard. For t > 5 Gyr, ah is typically 
10 — 30 AU. There is no strong pattern in eccentricities of 
ejected BH binaries - there are six with e ^ 0.8, seven with 
0.8 < e ^ 0.95, and only two with e > 0.95. The maxi- 
mum eccentricity of an ejected binary is e = 0.972. The BHs 
which are members of ejected binaries have a mean mass of 
10.98 Mq. This is more massive than the overall mean for 
BHs in Run 2, which have masses distributed uniformly in 
the range 8 ^ maH ^ 12 Mq. 

In addition to the 15 BH binaries, one triple BH system 
is ejected from Run 2, at r ~ 4100 Myr. This consists of a 
tight low-eccentricity binary {at = 8 AU, e = 0.376) with a 
single BH bound in a wider low-eccentricity orbit [ai, = 149 
AU, e = 0.370). 

Previous studies have demonstrated that binary BHs 
ejected from massive star clusters can have orbital prop- 
erties that would lead them to coalesce within a Hub- 
ble time due to the emission of gravitat ional radiation 
(see e.g.. iPortegies Zwart fc McMillaiil |200G| ). Such objects 
may therefore be possible candidates for detection by 
gravitational wave experiments. An approximate formula 
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for the time-scale for a BH binary to coalesce due to 
the emission of gravitational rad iation is given by (e.g., 
IPorteeies Zwart fc McMillanll200ol ): 



Tcoal ~ 3.2 X 



VmBH/ VAUy 



Gyr. 



(7) 



It is easy to show that none of the ejected BH binaries in 
Run 2 would merge within a meaningful time-scale (here 
we adopt ~ 12 Gyr, which is the approximate age of the 
Universe minus the delay of ~ 1.5 Gyr before the first BH 
binary ejections occur in Run 2). The most tightly bound 
ejected binary has Oj, = 4.2 AU and e = 0.839, while the 
most eccentric ejected binary has at = 8.2 AU and e = 
0.972. Even so, the orbital parameters of these objects are 
not vastly different from those which would lead to merging 
events on an interesting time-scale. For example, the binary 
with e = 0.972 would need a^, = 1.06 AU to merge in 12 
Gyr, while that with at = 4.2 AU would need e = 0.994. We 
consider this topic further in Sections 14.21 and [5] 

Recent large-scale A'^-body simulations have demon- 
strated comprehensively that when an intermediate- 
mass black hole (IMBH; mass ~ a few xIO^Mq) IS 
present in a massive star cluster, a central stellar den- 
sity and velocity cusp dev elops about this object (e.g., 
iBaumgardt. Makino fc Ebis uzaki 2004a!lb). It is natural to 
ask whether a similar cusp develops in Run 2, where a com- 
parable BH mass is concentrated in the cluster centre, but 
in the form of many relatively small objects rather than one 
massive object. 

Fig. [TT] summarises the structural and dynamical state 
of the stellar component of Run 2 at two output times; r = 5 
Gyr and r = 10 Gyr. These are late enough that any cusp 
should h ave had sufficient time t o form (see e.g., the time- 
scales in IBaumgardt et al.ll2004al lbh. The top panels in Fig. 
[TT]show the 3-dimensional radial mass density profile of the 
cluster at the two output times (solid circles). All luminous 
matter in the cluster was counted in each profile (i.e., BHs 
were excluded). The radial bins contain 50 stars for radii 
closest to the cluster centre, graduating to 100 stars, then 
500 and 1000 stars at increasingly large radii. 

For comparative purposes, we have also plotted depro- 
jected EFF models in these panels. These models are of the 
form of Eq. IA2I In calculating them, we used the values of 
/xo, a, and 7 observed from the projected brightness pro- 
file at the appropriate time, as described in Section O The 
maximum radial extent of the projected brightness profiles 
is marked in Fig. [Tl]by vertical dotted lines. Agreement be- 
tween the models and data is not necessarily expected be- 
yond these radii; in addition, tidal effects become important 
at the largest radii. For convenience (see below), we took 7 
in the deprojected models to be the closest integer value to 
that observed - that is, 7 = 4 at t = 5 Gyr, and 7 = 6 
at T = 10 Gyr. In each deprojected EFF model, the central 
surface luminosity density, /io, was converted to the volume 
luminosity density jo via Ea. lA2l For example, at r = 5 Gyr, 
we measured ^0 ~ 0.55magpc~^ — 51.05 L© pc~'^, which 
corresponds to jo ~ 4.59 Lq pc~^. To obtain a mass density 
from this value requires multiplication by a global mass-to- 
light ratio appropriate for the age and metal abundance of 
the cluster. We determined this empirically by fitting the 
deprojected EFF model to the measured data. The result- 
ing mass-to-light ratios [M/L = 1.33 at r = 5 Gyr, and 
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Figure 11. Radial mass density profiles and ID velocity disper- 
sion profiles for Run 2 (upper and lower panels, respectively) at 
cluster ages of 5 Gyr and 10 Gyr (left and right panels, respec- 
tively). In all panels, data from the simulation are marked with 
solid black circles. Only stars were used to derive these measure- 
ments - BHs were excluded. In the upper panels, the solid line 
indicates a deprojected EFF model fit using parameters set by 
those observed from the simulation at the appropriate output 
time. The vertical dotted lines indicate the maximum radius used 
in the construction of projected brightness profiles when obtain- 
ing these observations (see Section [Sjl. In each lower panel, the 
dashed line represents the velocity dispersion profile predicted by 
the EFF model plotted in the matching upper panel, while the 
solid lines represent the same models with rescaled central densi- 
ties to fit the data (see text). 



M/L — 2.01 at T = 10 Gyr) are a good match for those cal- 
culated by directly summing the mass and luminosity of all 
stars in the cluster, excluding BHs {M/L = 1.35 at r = 5 
Gyr, and M/L = 2.10 at r = 10 Gyr). We note that the 
assumption of a globally constant mass-to-light ratio is a 
reasonable approximation for Run 2 at these late times due 
to the relatively low degree of mass segregation amongst the 
stellar component of this cluster. 

It is clear from Fig. [TT] that no significant cusps are 
present in the cluster's stellar density profile at either time. 
At the most, it is possible that very marginal cusps exist, 
since the density profiles rise slightly above the EFF models 
(which have constant density cores) at the innermost few 
data points; however the significance of this "density excess" 
is very low. Certainly, striking cusps of the form of those 
observed by IBaumgardt et aL (|2004al lbh to develop about 
central IMBHs in clusters are not present. 

In the lower panels of Fig. 1111 we plot the ID stellar 
velocity dispersion as a function of radius at the two output 
times. The same radial bins as in the density profiles were 
used. Again, although the central regions of these profiles 
show some point-to-point scatter, there is no evidence at 
either time of a significant central velocity cusp analogous 
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to the type observed by iBaumgardt et al] (|2004al lbh when 
an IMBH is present. 

It seems hkely that the absence of a steUar density and 
velocity cusp about the central BH subsystem in Run 2 is 
due to the fact that scattered and ejected BHs are constantly 
moving through the region where a cusp would be expected 
to develop. This region is hence constantly being disturbed 
so that stars cannot settle into a stable distribution about 
the central concentrated mass as they can when only a single 
high-mass object is present. Such a process is similar to the 
destruction of cusps in galactic nuclei by supermassive black 
holes; except in that case a sin gle very massive binar y BH 
typically does the damage (e.g., Merritt &: CrudboOll ). 

Since there is no evidence from the 3D radial profiles 
for any large central density or velocity cusps, the projected 
profiles which it is possible to observe for real clusters will 
certainly show no evidence for any cusps. This is supported 
by the surface brightness profiles calculated at each output 
time in Run 2 to measure rc and 7, which exhibit constant 
density c ores as observed for the ma jority of LMC and SMC 
clusters (|Mackev fc Gilmore|[2003al lbh. 

Is there therefore another means by which we might in- 
fer observationally the presence of the significant central BH 
population in Run 2? The lower panels in Fig. 1111 together 
with Fig. [12] sketch the principles of one potentially viable 
method. In each of the lower panels in Fig. 1111 we plot the 
stellar velocity dispersion profile predicted by taking the pa- 
rameters (po,a,7) from the nicely-fitting deprojected EFF 
density model marked in the respective top panel. Since we 
chose integer values of 7, the predicted velocity dispersion 
profiles are analytic, and easily computed. That for 7 = 6 
is given by Eq. IA14I while the 7 = 4 case is the well known 
Plummer sphere: 



cr^(r) = 



2-KGpoa? 



(8) 



9a/1 + 7^ 



The resulting velocity dispersion profiles are plotted with 
dashed lines in the lower panels of Fig. 1111 Clearly they are 
a very poor fit to the measured profiles. However, simply 
rescaling the central density po so that the central veloc- 
ity dispersion predicted by the deprojected EFF model is 
consistent with the innermost measured data points results 
in rather nice fits (solid lines), at least out to large radii 
where the external tidal field begins to affect the stellar dy- 
namics. The required central densities at 5 and 10 Gyr are, 
respectively, « 3.4 and ~ 2.1 times those determined from 
the density profiles in the upper panels. This clearly implies 
that unseen matter (i.e., the BH population) is influencing 
the stellar dynamics in the cluster. By measuring the veloc- 
ities of stars in an extended cluster, we might therefore be 
able to infer the presence of a retained BH population. 

Such measurements are difficult. Apart from any tech- 
nical intricacies, we are limited to working in projection and, 
with present technology, to using line-of-sight velocities only. 
The left panel of Fig. [12] shows the dispersion in the line- 
of-sight velocities in Run 2 at r = 5 Gyr as a function of 
projected radius. As previously, all luminous matter in the 
cluster was counted in the profile, but BHs were excluded. 
The radial bins contain 50 stars for radii closest to the clus- 
ter centre, graduating to 100 stars, then 500 and 1000 star s 
at increasingly large projected radii. [Wilkinson eFal] (|2002l l 
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Figure 12. Dispersion in the line-of-sight velocities in Run 2 at 
T = 5 Gyr, as a function of projected radius. In the left panel, data 
from the simulation are marked with solid black circles. Only stars 
were used to derive these measurements - BHs were excluded. 
The dashed line represents the dispersion profile predicted by a 
Plummer sphere with central density taken from the profile in 
the upper left panel of Fig. 1111 The solid line represents the same 
model with a rescaled central density (see text). In the right panel 
dispersion profiles from two stellar sub-samples are plotted: solid 
dots are for a sample of 1000 upper main sequence stars, and 
crosses are for a sample of 150 red giant branch stars. The two 
models from the left-hand panel are also marked. 



provide an expression for crios(rp) in a Plummer sphere (their 
Eq. 27 and 28). In the simplest case of an isotropic mass- 
foUows-light model their expression reduces to: 



'^los('"p) 



TT Gpoa 



(9) 
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which has the same functional form as cr^{r), but a slightly 
smaller central value: crf^siO) ~ 0.88o"^(0). This suggests that 
the dynamical signature we observed from the lower panels 
of Fig. [11] should still be visible in projection, and indeed 
we find this to be the case. In the left panel of Fig. [12] we 
fit a model of the form Eq. [^to the data (solid line), again 
leaving po as a free parameter. In this case we find po — 
19.5Mq , which is very similar to the value required to fit the 
deprojected velocity data, but very different from the value 
implied from the radial density profile. We also computed 
a model using this latter value (po ~ 6.1Mq); this is the 
dashed line in the left panel of Fig. [T^] 

The difference between the two models is sufficiently 
large that it may be detectable in clusters using presently 
available technology. In the right panel of Fig. [T^] we plot 
(Jiasii^p) determined using two samples of stars randomly se- 
lected from Run 2. The first is a sample of 1000 upper main 
sequence stars grouped into six bins of 125 stars, while the 
second is a sample of 150 red giant branch stars grouped into 
five bins of 30 stars. Both samples clearly favour the model 
with the larger mass-to-light ratio. The red-giant sample 
is of a size which could feasibly be measured by a modern 
multi-object spectrograph such as VLT/FLAMES, although 
it must be borne in mind that the typical measurement er- 
rors in radial velocities observed with such a facility will 
be at least comparable in magnitude to the dispersion in a 
diffuse cluster (~ 1 — 2 kms~^) - a sophisticated analysis 
would be required to properly account for these. 

While somewhat crude, our results demonstrate that 
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in a cluster such as that modelled in Run 2, where there 
is a relatively large BH population present, the stellar dy- 
namics should imply the presence of significantly more mass 
than is evident from observations of the luminous compo- 
nent of the cluster. This arms us with a means of searching, 
albeit indirectly, for BH populations in massive LMC and 
SMC star clusters. Even so, we expect such observations to 
be extremely challenging due to the small velocity disper- 
sions involved, the necessity of working in projection, and 
the general sparsity (in terms of numbers of bright stars) of 
the extended clusters observed in the Magellanic Clouds. 

4.2 Runs 3 and 4: Strong mass segregation 

We next consider the pair of simulations labelled Run 3 and 
Run 4. These are both strongly primordially mass segregated 
clusters, created as described in Section 12.21 using a pre- 
evolution duration of Tms ~ 450 Myr. This duration was 
selected so that Runs 3 and 4 possess initial properties very 
similar to those observed for the very young, compact cluster 
R136 in the 30 Doradus complex in the LMC (see below). 

Like Runs 1 and 2, Runs 3 and 4 start with identical 
initial conditions, to the extent that they share the same 
random seed. Once again, the sole difference between them 
is that in Run 3 the natal BH kicks are set to be Ukick ~ 200 
kms"'^, whereas in Run 4 they are set to be zero - this 
results in /bh = and /bh = 1, respectively. 

The primary aim of Runs 3 and 4 is to try and follow 
the evolution of models which look more similar to the very 
youngest (r < 20 Myr) massive LMC and SMC clusters than 
do Runs 1 and 2. In particular, as discussed in Section [2. 2 1 a 
number of very young Magellanic Cloud clusters have been 
observed to be mass segregated to some degree. However, 
significant mass segregation does not develop in Runs 1 and 
2 until T ~ 100 Myr or so. In addition, the projected bright- 
ness profiles of Runs 1 and 2 (and particularly the structural 
parameters Tc and 7) do not resemble observed young LMC 
and SMC clusters until t ^ 20 Myr (e.g.. Fig. [5j. These 
differences mean that the observed early evolution of Runs 
1 and 2 may not accurately reflect the processes occurring 
in the youngest massive Magellanic Cloud clusters. 

Furthermore, in Run 2 we found that the BH population 
did not influence the structural evolution of the cluster until 
after the formation of the first BH binaries in the core at r ~ 
510 Myr. Since Fig. [T] shows that there is clearly evolution 
in the observed radius-age trend on time-scales shorter than 
this, it is important to examine whether it is possible for the 
BH population to become dynamically active earlier than 
seen in Run 2. One might naively expect this to occur if BHs 
are formed preferentially at the centre of a cluster, such as 
they would be in a primordially mass segregated object. 

It is important to first assess the suitability of the initial 
conditions we constructed for Runs 3 and 4 before moving 
on to an examination of the evolution of these Runs. One 
simple but useful indication is provided by the observed ini- 
tial structural parameters rc, 7, and po- The measured val- 
ues for Runs 3 and 4 are listed in Table [T] As described 
previously, these quantities are an excellent match for those 
determined for R136; see also Fig.|3] We note however, that 
R136 is nearly an order of magnitude more massive than our 
TV-body models. Scaling issues are discussed in Section [S] 

One of the major differences between Runs 3 and 4 and 
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Figure 13. Early evolution of the central density po for Run 1 
and Run 3. Run 1 has no primordial mass segregation while Run 
3 is strongly primordially segregated. Run 3 appears very similar 
to R136 at early times; however by a few tens of Myr its central 
density has decreased significantly and is a good match to LMC 
and SMC clusters of this age. Run 1 starts with a much lower 
central density, which it maintains throughout its early evolution. 
Together these two models span the observed density ranges for 
the youngest LMC and SMC objects. 



Runs 1 and 2 is that the former have very much larger cen- 
tral densities than the latter. This is simply due to the strong 
central concentration of the most massive stars in Runs 3 
and 4 as a result of the initial mass segregation. It is enlight- 
ening to examine the early evolution of the central densities 
in these differing models - this evolution is plotted for Runs 
1 and 3 in Fig. 1131 At the start of the simulation, the density 
of Run 3 is directly comparable to that of R136; however, 
as the early phase of severe mass loss due to stellar evo- 
lution begins, the central density drops rapidly so that by 
~ 10 Myr it matches the densities observed for other young 
LMC and SMC clusters. This rapid drop in central density 
implies that the central regions of the cluster expand dur- 
ing this early period of evolution, and indeed this is what is 
observed (see below). In comparison, the initial central den- 
sity of Run 1 is much lower than that of R136, and does not 
change much as the rapid early stellar evolution progresses. 
This is consistent with the fact that Run 1 shows little or no 
central expansion during this phase. Together, Runs 1 and 3 
span the range of central densities observed for the youngest 
LMC and SMC clusters - we are therefore confident of the 
applicability of our models in this regard. 

It is also possible to assess how well the primordial mass 
segregation generated in Runs 3 and 4 matches that ob- 
served in genuine young Magellanic Cloud clusters. We do 
this by performing simulated observations of the radial vari- 
ation of the stellar mass function (MF) in the models, and 
comparing the results t o those deterrn i ned from th e detailed 
observational studies of'Huntcr et al. ' ('1995', 'iggd) for R136 
(~ 3 Myr old); that of^le Grijs et aL (2002a) for NGC 1805 
(~ 10 Myr old) and NGC 1818 (~ 20 Myr old); and that of 
ISirianni et all (|2002h for NGC 330 (~ 30 Myr old). 

In performing the simulated observations, we follow the 
individual cluster studies as closely as possible. That is, we 
use the same projected radial bin widths and the same stel- 
lar detection limits within each such bin as were used in 
the original observational study. This ensures that our mea- 
surements are closely comparable to t hose obtained i n each . 
Consider, for example, the work of ISirianni et al.) l|2002h . 
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Figure 14. Mass and luminosity function slopes as a function of projected radius for various young massive LMC and SMC clusters, 
compared with results from simulated observations of TV-body Run 3. The plots have been reproduced to match those pres ented for each 
cluste r by th e original authors. Left: Mass function slope T as a function of projected radius in R136 in the LMC from lHunter et al.l 
lll995l . ll99^ . compared with Run 3 at age 3 Myr. Centre: Luminosity function slopes /3 as a function of projected radius for NGC 
1805 and NGC 1818 in the LMC from lde Griis et all l l2002al V compared w i th Ru n 3 at age 15 Myr. Right: Mass function slope V as 
a function of projected radius in NGC 330 in the SMC from lSirianni et al. 1 liooi), compared with Run 3 at age 30 Myr. 



These authors used five annuli of 5" width to span the range 
— 25" in projected radius in their study of NGC 330, fol- 
lowed by ten annuli of 10" width to span the range 25—125" 
in projected radius. Ultimately, however, they decided to 
employ a maximum projected radius of 95" for their MF 
calculations, due to a significant contaminating population 
of field stars. Within their radial bins, the stellar mass lim- 
its used to calculate the MF were from the top of the main 
sequence in this cluster (~ 7M0) to the 50 per cent com- 
pleteness level - at ~ 1.3 M© in the cluster centre, decreas- 
ing to ~ 0.8 Mq at a projected radius of 25" and beyond. 
When measu ring our model cluste r for a comparison with 
the results of ISirianni et all (|2002h . we took the data from 
the output time nearest to 30 Myr, projected the positions 
of all stars onto a plane, converted the projected radial scale 
from parsecs to arcseconds by applyin g the SM C distance 
modulus of 18.85 assumed by. Sirianni et al. (2002;), and then 
a pplied exactly the sam e bin widths and mass limits per bin 
as lSirianni et al. I(l200t ) to o btain the stellar samples for MF 
fitting. ISirianni et al. rnooi) corrected their star counts for 
completeness variations, so we assumed 100 per cent com- 
pleteness in each radial bin. At our chosen output time, the 
core radius and central density of our model are within ~ 15 
per cent of the values measured for NGC 330 (see Figs. 1131 
and I15p . so the bins are sampling equivalent regions in the 
cluster. 

The results of our simulated observations may be seen 
in Fig. 1141 along with the original measurements for R136, 
NGC 1805 and 1818, and NGC 330. For R136 and NGC 
330, the mass functions are represented by C,{m), which is 
the number of stars per logarithmic mass interval, as op- 
posed to the mass function ^(m) defined in Section [2.11 If 
the mass function f (m) has a power- law exponent — a^, then 
the mass function C,{Tn) will have exponent Ti — —ai + 1. 
In the case of NGC 1805 and 1818, we calculate and fit lu- 
minosity fun ctions (LFs) rather th an MFs, so as to match 
the work of Ide Griis et al] (|2002al ). The assumed form of 
the LFs, defined here as the number of stars per logarithmic 
luminosity bin, is a power law of slope f}. 



From Fig.[T4]it is clear that the results obtained via our 
simulated observations of Run 3 are generally an excellent 
match for those measured for the real LMC and SMC clus- 
ters. The greatest diflterences occur for the innermost radial 
bin of R136, and the outermost radial bins of NGC 330. The 
former discrepancy suggests that the very centre of Run 3 
(within ~ 0.1 pc) may initially be somewhat more mass seg- 
regated than R136, although we note that the rest of our 
measurements are highly consistent with the real observa- 
tions of R136, and that the region taken by iHunter et al.l 
(| 19961 ) to obtain their measurement at the innermost radius 
is extremely crowded with bright stars. The latter discrep- 
ancy may be related to the necessity for significant field 
st ar decontamina t ion in the outermost regions of NGC 330 
by ISirianni et al.l (|2002l ) - again, we note that the majority 
our measurements of Run 3 are in excellent agreement with 
those obtained by these authors for NGC 330. 

Overall, these results are strongly suggestive that the 
initial conditions we constructed for Runs 3 and 4, and in 
particular of the algorithm we developed to generate the 
primordial mass segregation in these models, are valid. We 
note however, that we are not able to place any similar ob- 
servational constraints on the initial velocity distributions 
in these models. As an aside, it is extremely interesting to 
observe the progression of the radial mass/luminosity func- 
tion profile of Run 3 from age t — 3 Myr to 30 Myr, in 
comparison with the profiles observed for four genuine LMC 
and SMC clusters. While there has previously been noth- 
ing to link the measurements of these four objects, the early 
evolution of Run 3 clearly demonstrates that a cluster which 
initially possesses a core radius, central density and radial 
MF profile very similar to that of R136 can evolve to look 
very much hke NGC 1805 and NGC 1818 after 15 Myr and 
then further to look like NGC 330 after another 15 Myr, 
simply via internal cluster dynamical processes under the 
influence of rapid mass loss due to stellar evolution. 

In Fig. [T51 we show the evolution of Runs 3 and 4 across 
the radius-age and 7-age planes. Unlike Runs 1 and 2, both 
Runs 3 and 4 exhibit dramatic core expansion right from the 
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Figure 15. Structural evolution of Af-body Runs 3 and 4. Both possess significant primordial mass segregation; the only difference 
between them is the BH retention fraction (/bh = and 1, respectively). Left panel: Evolution of rc, observed as described in Section 
|3] Both models experience significant expansion over the first ~ 200 Myr of evolution, due to the early phase of severe mass loss due to 
stellar evolution. This is in contrast to Runs 1 and 2, where the mass loss was spread throughout the cluster rather than being centrally 
concentrated. Subsequently, Run 3 begins to relax dynamically and slowly contracts, whereas the BH population retained in Run 4 
becomes dynamically active, leading to further core expansion in this model. By Tmax = 10 Gyr, the core radius for Run 4 has moved 
off the top of the plot, to Tc ~ 11 pc. Right panel: Evolution of the power-law fall-off, 7, again observed as described in Section [3] 
Both models develop increasingly steep 7 values as they evolve; however, that for Run 3 reaches a plateau once the core expansion in 
this model ceases. In contrast. Run 4 develops a very steep fall-off in its outer regions. 



beginning of their evolution. This is in response to the early 
phase of severe mass-loss due to stellar evolution, which in 
Runs 3 and 4 is highly centrally concentrated because of 
the primordial mass segregation. Fig. [16] shows the radial 
distributions of all supernovae in Runs 1 and 3 - the more 
centralized location of these events in Run 3 compared with 
Run 1 is clearly evident. The central concentration of the 
mass-loss, together with the high central density in Runs 
3 and 4, means that the amount of heating per unit mass 
lost is maximised in these models, hence leading to the ob- 
served dramatic core expansion. This core expansion is at 
least partly responsible for the rapid decrease in the central 
density of Run 3 which we noted in Fig. [13] (the demise of the 
most massive cluster stars also contributes to this decrease) . 

It is interesting that Runs 3 and 4 do not undergo 
an early core collapse despite their high central densities. 
Early core collapse in a massive cluster may lead to a run- 
away merger event, which is one possible formation chan- 
nel for a central intermediate - mass b l ack hole (IMBH) (e.g., 
Portcgics Zw art fc McMiUanI I2OO2I : jPortegies Zwart et all 
2004). Previous work has demonstrated that clusters with 
a very short initia l median relaxation time are suscep tible 
to early collapse - jPorteeies Zwart fc McMillan! (|2002l ) sug- 
gest trh < 25 Myr. It is not clear whether a similar thresh- 
old is applicable to our primordially mass segregated models. 
These have very much longer initial median relaxation times 
(trh ~ 1.2 Gyr), but rather short central relaxation times 
(tro ~ 9 Myr). It is possible that expansion of the clus- 
ter core due to mass-loss from rapid stellar evolution acts 
against the tendency of the core to collapse more strongly 



in our models than in previous models, due to the initial 
preferentially central location of many massive stars. 

By T « 100 Myr the rate of mass-loss from stellar evo- 
lution has significantly decreased, and by r ~ 200 Myr the 
core expansion in Runs 3 and 4 has largely stalled. Even 
though both models initially contain identical stellar popu- 
lations. Run 3 loses more mass up to this point than does 
Run 4, because the 198 BHs in Run 3 escape from the clus- 
ter immediately after formation, whereas those in Run 4 are 
retained. This difference is reflected in the larger degree of 
early core expansion observed in Run 3 compared to Run 4. 
Up until an age of roughly a few hundred Myr, Run 3 closely 
traces the observed upper envelope of the radius-age trend. 

In both models, the early mass-loss and core expan- 
sion is accompanied by a significant steepening in the outer 
power-law fall-off. This is again in contrast to the evolution 
observed for Runs 1 and 2 during the early mass-loss phase, 
where 7 remains essentially constant. Similarly to the core- 
radius evolution in Runs 3 and 4, the steepening of 7 stalls 
beyond r « 100 Myr in these models, once the rate of mass- 
loss has decreased. Furthermore, the evolution of 7 up to 
this point is slightly different in the two models, due to the 
retention of BHs in Run 4 and their expulsion in Run 3. 

Thereafter, Runs 3 and 4 begin to evolve quite differ- 
ently. Run 3 progresses in exactly the fashion of Run 1 
- the cluster settles into a quasi-equilibrium state where 
the dynamical evolution is dominated by two-body relax- 
ation processes, leading to the gradual development of mass 
stratiflcation. Because Run 3 is far less dense than Run 1 
by this stage, its relaxation time-scale is much longer. By 
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Figure 16. Radial distributions of supernova explosions in Run 
1 (solid line) and Run 3 (dashed line). All supernova explosions 
in both models have occurred by r 50 Myr. The significantly 
greater central concentration of the mass loss in Run 3 compared 
to Run 1 is quite evident. 



Tmax = 10.27 Gyr, only 3.1 integrated median relaxation 
times have elapsed in this model, compared with 4.7 median 
relaxation times in Run 1 at the same age. Hence, while Run 

3 is evolving towards core collapse when the simulation is 
terminated, we would expect it to enter this phase at a much 
older age than observed for Run 1. 

In contrast, at r ~ 750 Myr, core expansion resumes in 
Run 4. This continues until the end of the simulation, which 
is terminated at Tmax = 10.0 Gyr. By this time, the core 
radius of Run 4 has evolved off the top of Fig. 1151 to reach 
almost 11 pc. This is roughly comparable in size to the core 
radii observed for the most extended old Magellanic Cloud 
clusters, such as Reticulum in the LMC and Lindsay 1 in 
the SMC (Mackey et al. 2007, in prep.). 

The second, prolonged, period of core expansion in Run 

4 is due to the dynamical activity of its retained BH pop- 
ulation. The evolution of this BH subsystem, illustrated in 
Fig. 1171 is qualitatively identical to that which we observed 
in Run 2. The BHs, once formed, sink via dynamical friction 
and begin to accumulate at the centre of the cluster (panels 
b and c) . The density of this central BH subsystem increases 
until it becomes sufficiently high as to initiate the creation 
of stable BH binaries in three-body interactions (panel a). 
The first such object is formed in Run 4 at ~ 570 Myr. Sub- 
sequently, these binaries undergo superelastic collisions with 
other BHs in the cluster centre, which leads to the harden- 
ing of the binaries, the scattering of BHs beyond r^, and the 
ejection of BHs (panels d and e). These processes in turn 
result in the observed long-term core expansion. 

One intriguing aspect of the evolution of the BH subsys- 
tem in Run 4 is that even though this cluster was strongly 
primordially mass segregated, so that the majority of the 
BHs were created in its very inner regions (cf. Fig. I16|) . the 
first BH binary does not form until a similar time as that 
in the non- segregated Run 2. From Fig. 1171 it is also clear 
that the peak central BH density occurs at a similar time 
as it does in Run 2. Fig. [T5] shows an expanded view of 
the early evolution of the BH Lagrangian radii in Run 4, 
with those for Run 2 plotted for comparison. As expected, 
the majority of BHs are formed near the centre of the clus- 
ter - the BH Lagrangian radii are initially much smaller 
than in Run 2. However, unlike Run 2, Run 4 suffers signif- 
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Figure 17. Properties of the BH population in Run 4 as a func- 
tion of time: (a) the number of single BHs (upper line) and binary 
BHs (lower line) in the cluster; (b) the number of BHs within 
the shells encompassed by the stellar Lagrangian radii r ^ Ri% , 
^1% < ^ ^ ^5% I ^iid > ^10% (ths uppermost to lowermost 
lines, respectively, at the right of the plot); (c) the black hole 10%, 
25%, 50% and 75% Lagrangian radii (respectively, the innermost 
to outermost lines) ; (d) the cumulative numbers of escaped single 
BHs (upper line) and binary BHs (lower line), along with fits of 
the form Ne = Aq + Ait — j42''"logT (dashed lines) and the cu- 
mulative numbers of escaped single and binary BHs measured for 
Run 2 (dotted lines); and (e) the radial positions of three typical 
BHs. The vertical dotted line indicates r = 750 Myr, the approx- 
imate time when core expansion due to BH activity begins. The 
evolution of rc is plotted (dot-dashed line) in panel (e). Note the 
different axis scales on either side of panels (a) and (d). 



icant early core expansion due to the rapid stellar evolution 
phase. The BH subsystem does not escape this - the early 
centrally concentrated mass-loss is severe enough that the 
resulting expansion overcomes the natural tendency of the 
BHs to sink to the cluster core. This is reflected in the swift 
outward movement of the BH Lagrangian radii in Run 4, un- 
til the mass-loss phase is mostly complete around r ~ 100 
Myr. Subsequently, the BHs do begin to sink to the cluster 
centre via dynamical friction, and the evolution of the BH 
Lagrangian radii in Run 4 closely follows that in Run 2. 
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Figure 18. Evolution of the 10%, 25% and 50% BH Lagrangian 
radii in Run 4 (solid lines). The evolution of the same radii in Run 
2 are plotted for comparison (dotted lines). This plot clearly shows 
that the BH subsystem in Run 4 expands at early times along with 
the rest of the core, in response to the rapid mass-loss from stellar 
evolution which is occurring. Once this phase is mostly complete 
(t « 100 Myr), the evolution of the BH subsystem is very similar 
to that in Run 2. 



Figure 19. Separations and eccentricities of the ejected BH bi- 
naries in Run 4 as a function of cluster age. Eccentricity is repre- 
sented by point style: BH binaries with e ^ 0.8 are asterisks, those 
with 0.8 < e ^ 0.95 are open circles, while those with e > 0.95 are 
filled circles. The ejected BH binaries from Run 2 are plotted for 
comparison (small crosses). The arrow marks the ejection time of 
one additional Run 4 BH binary, which has separation ai, = 56.4 
AU and eccentricity e = 0.609. 



This result suggests that, contrary to expectations, the 
BH population in a primordially mass segregated or cen- 
trally concentrated cluster does not become dynamically 
active at significantly earlier times than does an identical 
population in a non-segregated cluster. In turn, this implies 
that the evolution in the radius-age trend observed on time- 
scales shorter than « 500 Myr is not due to the influence of 
a BH population, unless such populations are comprised of 
BHs with masses significantly in excess of 10 Mq. Instead, 
the early evolution of the radius-age trend most probably 
reflects centrally concentrated mass-loss in dense clusters 
due to rapid stellar evolution. Never the less. Runs 3 and 
4 clearly demonstrate that this process cannot alone propel 
clusters to the upper right corner of Fig. [T] since it oper- 
ates on a time-scale which is much too short. Our A'^-body 
models which possess core radii larger than ~ 6 pc after a 
Hubble time of evolution do so only because they have ex- 
perienced prolonged core expansion due to the activity of a 
retained BHs, irrespective of whether they also experienced 
core expansion at ages r < 100 Myr. 

As with Run 2, by Tmax = 10.0 Gyr there is still a 
significant BH population remaining in Run 4: 95 single BHs 
and 2 binary BHs. In fact, this population is rather larger 
than that in Run 2 at the same age. From Fig. I17t l. it is clear 
that, while the cumulative numbers of escaped single and 
binary BHs in Run 4 follow the same functional dependence 
on age as in Run 2, they are, at all times, smaller than those 
in Run 2. That is, the rates of escape of BHs are always 
somewhat lower in Run 4 than they are in Run 2. 

We attribute this variation to the different overall struc- 
tures of Runs 2 and 4 when their respective BH populations 
become dynamically active. In Run 2 the core radius and 
central density remain nearly constant from the beginning 
of the simulation until this point (cf. Figs. [S] and I13p : in 
contrast, Run 4 undergoes significant early core expansion. 
By T = 500 Myr, Run 4 is a considerably more diffuse clus- 
ter than is Run 2. The shallower gravitational potential in 
Run 4 affects the distribution of the BHs within this cluster 
(cf. Fig. I17fc and Fig. [T];). This leads to a slower interaction 



rate between BHs in Run 4 than in Run 2, and hence the 
observed lower BH escape rates. The same effect is primar- 
ily responsible for the BH escape rate in a model cluster 
decreasing as the core radius increases (see Section |4Tj, al- 
though in that case the decreasing size of the BH population 
contributes in addition. 

The more diffuse nature of Run 4 also affects the prop- 
erties of the ejected BH binaries in this simulation com- 
pared with those in Run 2. In Run 4, there are 12 BH bi- 
naries ejected over the course of the simulation. Their sepa- 
rations and eccentricities are plotted in Fig. 1191 along with 
the ejected BH binaries from Run 2 for comparison. Because 
Run 4 is always more diffuse than Run 2 at times when bi- 
nary BHs exist, these objects are more easily ejected in Run 
4 than in Run 2. Hence, the ejected BH binaries in Run 4 are 
generally not as tightly bound as those in Run 2. This can 
be seen in terms of the binary separations in Fig. 1191 which 
are typically larger for the ejected binaries in Run 4 than 
for those in Run 2 at similar times. In addition, the ejected 
BH binaries in Run 4 are typically less eccentric than those 
in Run 2 - of the 12 ejected Run 4 BH binaries, only one 
has e > 0.95, while there are three with 0.8 < e ^ 0.95 and 
eight with e ^ 0.8. The maximum eccentricity of an ejected 
binary is e = 0.981, while the minimum is e = 0.225. As in 
Run 2, the members of ejected binaries are typically more 
massive than the average mass for the full BH population - 
the mean mass of members in escaping binaries in Run 4 is 
11.20 Mq . In Section |4]T] we showed that none of the ejected 
BH binaries from Run 2 would merge due to the emission of 
gravitational radiation within a Hubble time (Eq. [7J . This 
is also, unsurprisingly, the case in Run 4. 

Fig. I20l shows the evolution of the ratio of core radius to 
projected half-light radius, rc/rh,i, for Runs 3 and 4. Runs 1 
and 2 are also plotted, for comparative purposes. The initial 
value of rc/rh,i for Runs 3 and 4 is significantly smaller than 
that for Runs 1 and 2, reflecting the primordial mass segre- 
gation in these models. However, the early core expansion 
in Runs 3 and 4 results in a rapid and significant increase in 
rc/rh,i- Overall, Runs 3 and 4 better match the observed dis- 
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Figure 20. Evolution of the ratio of core radius to projected 
half-liglit radius r-c/r^ ; for A^-body Runs 3 and 4, compared with 
measurements for LMC and SMC clusters. The evolution of Runs 
1 and 2 is also plotted, for comparative purposes. 



tribution of young and intermediate-age Magellanic Cloud 
clusters than do Runs 1 and 2. By the end of Run 3, rc/rh,i is 
beginning to decrease as two-body relaxation begins to dom- 
inate in this model; however, for the majority of this Run 
fcl'<'K,i ~ 0.7. This demonstrates that a large observed value 
of rc/rh,i in a physically old star cluster need not reflect the 
presence of a central massive body (such as an IMBH) or 
a central accumulation of many less-massive bodies (such 
as stellar-mass BIfs), but rather may reflect the fact that 
such a cluster is not very dynamically old. Run 4, which un- 
dergoes prolonged core expansion due to its BH population, 
finishes with rc/rh.i ~ 0.8, matching Run 2 closely. 

4.3 Runs 4a and 4b: Variable mass segregation 

As described at the beginning of Section UJ our four primary 
simulations cover the extremes of the parameter space we are 
investigating (spanned by < /bh 1 and < Tms < 450 
Myr), and are therefore expected to represent the extremes 
of cluster evolution induced by variation of these particu- 
lar initial conditions. However, it is important to sample 
intermediate values of both /bh and the degree of primor- 
dial mass segregation to check that the parameter space is 
well behaved and that the models evolve as we expect (i.e., 
intermediate between the extremes of Runs 1-4). 

With this in mind, we completed two additional simula- 
tions with /bh = 1, and degrees of primordial mass segrega- 
tion spaced between that for Run 2 and that for Run 4. Be- 
cause the initial conditions for these new models were taken 
from two different output times during the pre-evolution of 
Run 4, at Tms = 115 Myr and 225 Myr, we denote them as 
Runs 4a and 4b, respectively. We only ran these models to 
Tmax ~ 4.3 Gyr, as this was more than sufficient to observe 
the progress of the clusters' evolution. 

The initial properties of Runs 4a and 4b are listed in 
Table [T] As expected, the values of rc, 7, and po all lie be- 
tween those of Runs 2 and 4. The longer the duration of the 
pre-evolution, the smaller the initial value of rc, the higher 
the initial value of po, and the flatter the initial value of 
7. This latter, in particular, is expected due to the devel- 
opment of a cor e-halo structu re in a dynamically evolving 
cluster (see e.g.. ISpitzerlll987h . 

The core radius evolution of Runs 4a and 4b is illus- 
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Figure 21. Core radius evolution of Af-body Runs 4a and 4b, 
with the evolution of Runs 2 and 4 plotted for comparison. Runs 
4a and 4b, with pre-evolution durations of Tms = 115 Myr and 
Tms = 225 Myr, possess primordial mass segregation of degrees 
intermediate between those of Runs 2 and 4. This initial condition 
is the only difference between each of the four plotted Runs - all 
four form 198 BHs and have a BH retention fraction of /bh = 1. 



trated in Fig. 1211 Just as with Run 4, these two models 
undergo two main stages of core expansion. The first, last- 
ing until a little after ~ 100 Myr, is in response to the early 
rapid stellar evolution. The second, which begins around 
r ~ 600 — 800 Myr is due to the influence of the retained 
BH population. In between these two phases, there is a pe- 
riod during which the core radius is roughly constant. 

As expected, the core radius evolution seen for both 
Runs 4a and 4b lies between the limits defined by Runs 

2 and 4. The amount of early expansion apparently varies 
directly with the degree of primordial mass segregation - the 
more mass segregated a cluster, the larger the core expansion 
seen at ages less than ~ 100 Myr. From Fig. 1151 for Runs 

3 and 4, we saw that the amount of mass lost during the 
early period of rapid stellar evolution also influences to some 
extent the degree of the observed core expansion. However, 
all four models in the present case were speciflcally designed 
to possess identical populations of massive stars and retained 
BHs, and all therefore lose essentially identical amounts of 
mass due to stellar evolution at early times. The variation 
in the core expansion observed during this phase in Fig. 
[21] therefore cannot be caused by differing amounts of mass- 
loss and must solely be a result of the different initial cluster 
structures. More centrally concentrated mass- loss, in a more 
tightly-bound cluster core, clearly leads to a greater degree 
of core expansion during the early period of a cluster's life. 

After the first stage of core expansion is complete in the 
four model clusters, their core radii remain roughly constant 
until the BH populations have accumulated at the cluster 
centres and started to form BH binaries, after which point 
the second phase of core expansion begins. From Fig. 1211 
the rates of observed expansion in this second phase are 
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Figure 22. Evolution of the 10%, 25% and 50% BH Lagrangian 
radii in Runs 4a (upper panel, solid lines) and 4b (lower panel, 
solid lines). The evolution of the same radii in Run 4 are plotted 
for comparison in both panels (dotted lines). The BH subsystems 
in Run 4a and 4b both expand at early times in response to the 
rapid mass-loss from stellar evolution which is occurring, although 
the expansion is greater in the more heavily mass segregated Run 
4b. Once the early mass-loss phase is mostly complete (r ^ 100 
Myr), the evolution of the BH subsystems are very similar to that 
in Run 4 (and Run 2 - cf. Fig.lTSll. 



approximately equivalent in all four clusters, so that the 
tracks on the radius-age plane run close to parallel for the 
remainder of the simulations. 

Fig. [21] shows that the second stage of core expansion 
begins at an approximately equivalent time in each of the 
four models. We already noted this fact for Runs 2 and 4 in 
Section 14.21 and concluded that in a primordially mass seg- 
regated cluster the BH population does not become dynam- 
ically active significantly earlier than in a non-segregated 
cluster, because the strong expansion experienced by the 
mass segregated cluster at early times affects the BH pop- 
ulation sufficiently strongly to negate the natural tendency 
of the BHs to sink to the cluster centre. In Fig. [22] we plot 
the evolution of the BH Lagrangian radii in Runs 4a and 
4b, with those for Run 4 plotted for comparison. It is clear 
from this plot that even though the BHs tend to form closer 
to the cluster centres in more primordially mass segregated 
models, these models also suffer greater degrees of early ex- 
pansion, hence delaying the central accumulation of BHs. 
This results in very similar evolution of the BH Lagrangian 
radii in Runs 2, 4, 4a, and 4b after the early rapid stel- 
lar evolution phase is complete, and leads to the formation 
of the first BH binaries at very similar ages - 510 Myr, 570 
Myr, 540 Myr, and 460 Myr, in the four models respectively. 
Given that this is by nature a stochastic process, the agree- 
ment between these times for four models with such a wide 
range of early structural evolution is quite close. This ob- 
servation strengthens our earlier conclusion that primordial 



mass segregation in a cluster does not lead to the earlier de- 
velopment of a dynamically active BH subsystem compared 
to a cluster which is not primordially mass segregated. 



4.4 Run 5: Intermediate BH retention 

Together with Runs 4a and 4b, we computed one additional 
model possessing properties intermediate between those of 
our four primary Runs. In this case, instead of intermediate 
degrees of primordial mass segregation, we set up the sim- 
ulation (labelled Run 5) so that /bh ~ 0.5. Its initial con- 
ditions were otherwise identical to to those in Runs 3 and 4 
(i.e., strong primordial mass segregation set by Tms = 450 
Myr). One aim of Run 5 is to check that, as should be ex- 
pected, its core radius evolution lies between that observed 
for Run 3 (where /bh = 0) and that observed for Run 4 
(where /bh = 1). More importantly however, this model ex- 
plores whether the extreme case that /bh ~ 1 is necessary 
for significant core expansion to occur, or if such expansion 
can still develop in a system which loses a large fraction of 
its BHs at formation. We set the duration of Run 5 to be 
roughly a Hubble time (rmax = 10.06 Gyr) so that we could 
observe the full long-term core evolution of this model. 

To generate a retention fraction of roughly 50 per cent 
in Run 5, we examined the formation of BHs in Run 4 with 
the aim of determining a suitable distribution of natal kicks. 
More specifically, we calculated the potential energy ()7bh) 
of each given BH at the moment of its form ation in Run 4 , 
and estimated the escape velocity v^sc = -y^— 2(7bh / ^bh. 
Under the assumption that the inherited kinetic energy of 
the BH at formation (i^BH) does not contribute to its ejec- 
tion, this escape velocity corresponds to the minimum natal 
kick required to remove the BH from the cluster. However, 
this assumption is not always justified - for example, the 
natal kick may be in the same direction as the inherited 
motion of the BH, in which case the minimum required ve- 
locity would be significantly smaller than the origin al esti- 
mate, and closer to v^sc = V-2(f^BH + A'bh) / mBH. 

The two calculated distributions of BH escape velocities 
in Run 4 may be seen in Fig. 1231 The upper left panel is the 
binned distribution neglecting the inherited kinetic energy 
of the BH, while the upper right panel is the binned dis- 
tribution under the assumption that the kick velocity is in 
the same direction as the inherited velocity. The lower panel 
shows the two distributions in cumulative form. The most 
tightly bound BHs require natal kicks of order ~ 12 kms~^ 
to escape the cluster, while the least tightly bound BHs re- 
quire only tiny natal kicks to escape. Note from the upper 
right panel that two BHs are unbound from the cluster at 
their formation - this is a result of their progenitor stars 
becoming unbound shortly before exploding as supernovae, 
because of the recent rapid mass-loss occurring in the pro- 
genitor star as well as the violent fiuctuations occurring in 
the cluster potential due to mass-loss from other stars. 

For interest's sake, we also calculated the same distri- 
butions for Run 2, which, in contrast to Run 4, was not pri- 
mordially mass segregated. The distributions for this model 
are plotted in the lower panel of Fig. [23] as dashed lines. As 
might be expected, BHs formed in the mass segregated Run 
4 are significantly more tightly bound than are BHs formed 
in the non-segregated Run 2. Hence, the initial structure of 



Core expansion in massive star clusters 25 



0.1 



PE only 




PE + KE 




1 

S 0.8 - 




0.6 



0.4 



, , , 1 ,ni, ■■.■■^^^H 
5 10 5 


10 






- N-body Run 4 / 

! / 
f J 






/ /j 












. . 1 r— r--^ 1 1 1 1 1 1 1 







5 10 

Figure 23. Calculated distributions of escape velocities at forma- 
tion for all 198 BHs in Run 4. The upper left panel is the binned 
distribution assuming Vesc = rriBH, while the upper 

right panel is the binned distribution if the inherited kinetic en- 
ergy is also included so that Uesc 

The distributions may be interpreted as the minimum natal kicks 
required to remove the respective BHs from Run 4 under the two 
different assumptions listed above. The lower panel shows the two 
distributions in cumulative form (solid lines, where the distribu- 
tion including i^BH is to the left of that where ii'BH is neglected) . 
Also marked are the equivalent distributions for Run 2, which has 
no primordial mass segregation (dashed lines). 



a cluster can have some efltect on the retention fraction of 
BHs. We discuss this issue in more detail in Section (5] 

The distributions observed for Run 4 in Fig. [23] allowed 
us to determine a suitable distribution of natal kicks in Run 
5 in order to set /bh ~ 0.5. We did this by noting that the 
retention fraction /bh is approximately the mean probabil- 
ity of retention calculated over the full BH subsystem - that 
is, /bh ~ X/ (retain) / Nbh where A^bh is the number of 
BHs in the subsystem and Pi (retain) is the probability that 
the i-th BH will not be ejected by the natal kick it receives 
at formation. For simplicity, we set the natal kicks in Run 5 
to be selected randomly from a uniform distribution spread 
between Wkick = kms~^ and Wkick = tik.max kms~^. In this 
case, for the i-th BH at formation the retention probability 
is given by Pi(retain) = P(ukick,i < t'osc.i) = Ucsc,i / ^ik.max 
if ^ik.max > Vcsc,i, OT Unity Otherwise. Assuming that iik.max 
is larger than Wcsc,i for all BHs in the subsystem under con- 
sideration, in order to obtain a given retention fraction we 
require Uk,max = X] ^<=s<:.i / (/bhNbh)- For Run 5 we have 
that A^BH = 198, and require that /bh ~ 0.5, and we com- 
pute the sum using the distributions displayed in Fig. 1231 
We found that in this scenario Uk^max ~ 17.5 kms~^, de- 
termined by adopting the mean of the results for the two 
measured distributions (i.e., with and without the inherited 
BH kinetic energy). No physical meaning should be read into 
our choice of a uniform kick distribution - we selected it here 
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Figure 24. Core radius evolution of Af-body Run 5, with the 
evolution of Runs 3 and 4 plotted for comparison. All three models 
possess identical initial conditions, including strong primordial 
mass segregation (Tms = 450 Myr) and the formation of 198 BHs. 
The only difference between them is the BH retention fraction, 
which is zero for Run 3, unity for Run 4, and approximately 0.485 
for Run 5 — i.e., in Run 5 96 BHs are still present after 100 Myr. 



simply for convenience. The above process could easily be 
generalized to any desired probability distribution. 

With the natal kick distribution described above im- 
plemented in Run 5, as expected we observed a significant 
number of BHs escaping from the cluster shortly after their 
formation. All 198 BHs in the simulation are created by 
r = 10 Myr, and the first BH escapes occur at r = 14 
Myr. By r = 100 Myr, 102 BHs have left the cluster, and 
the escape rate has dropped to approximately one BH ev- 
ery 100 Myr. Subsequently, BHs appear to be leaving the 
cluster due to straightforward relaxation and evaporation 
processes - by the time of the creation of the first stable 
BH binary at r « 1350 Myr, a further 9 BHs have been 
removed. Hence, we estimate the BH retention fraction in 
this realization of Run 5 to be /bh ~ 0.485, but it could be 
as low as /bh ~ 0.440 depending on whether BHs escaping 
between r ~ 100 — 1350 Myr are included in the definition. 
We note that re-running this simulation with a new random 
seed would result in a different retention fraction, since, in 
contrast to all our previous models, the ejection of BHs due 
to natal kicks is a stochastic process in Run 5. 

The evolution of Run 5 on the radius-age plane is plot- 
ted in Fig. 1241 along with the evolution of Runs 3 and 4 for 
comparison. All three models possess identical initial condi- 
tions - the only difference between them is the BH retention 
fraction. At all times, the core radius of Run 5 lies in be- 
tween those measured for Runs 3 and 4. During the first, 
early, phase of core expansion, this is a consequence of the 
intermediate BH retention fraction in Run 5 - this model 
loses less mass than Run 3 but more than Run 4. The sec- 
ond phase of core expansion, due to BH activity, begins at 
T ~ 1400 Myr in Run 5. This is noticeably later than in 
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Figure 25. Early evolution of the 10%, 25% and 50% BH La- 
grangian radii in Run 5 (solid lines). The evolution of the same 
radii in Run 4 are plotted for comparison (dotted lines). The BH 
subsystem in Run 5 expands significantly at early times primarily 
due to the non-zero natal kicks, although the BHs do also share 
in the general expansion of the cluster due to the rapid mass-loss 
from stellar evolution which is occurring during this period. 
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Figure 26. Properties of the BH population in Run 5 as a func- 
tion of time: (a) the number of single BHs (upper line) and binary 
BHs (lower line) in the cluster; (b) the cumulative numbers of es- 
caped single BHs (upper line) and binary BHs (lower line); and 
(c) the black hole 10%, 25%, 50% and 75% Lagrangian radii (re- 
spectively, the innermost to outermost lines). The vertical dotted 
line indicates r = 1400 Myr, the approximate time when core ex- 
pansion due to BH activity begins. Note the different axis scales 
on either side of panels (a) and (b). 



the previous models; furthermore, the rate of expansion is 
clearly not as great as is observed in Run 4 where twice as 
many BHs are retained. 

The evolution of the BH population in Run 5 is illus- 
trated in Figs. [25] and 1261 The first of these shows the very 
early evolution of the 10%, 25% and 50% BH Lagrangian 
radii, compared to Run 4. Upon formation of the BHs, the 
BH Lagrangian radii are initially identical in Runs 4 and 



5 because they share identical initial conditions and ran- 
dom seeds. However, the Lagrangian radii in Run 5 immedi- 
ately undergo a large degree of expansion. This is primarily 
in response to the non-zero natal kicks given to the BHs, 
although the BH subsystem does share in the general ex- 
pansion of the cluster which is occurring during this period 
due to the rapid mass-loss from stellar evolution. For com- 
parison, the Run 4 Lagrangian radii are expanding only in 
response to this mass-loss. The Lagrangian radii in Run 5 
exhibit large bumps at early times - these are most promi- 
nent in the 25% radius at r < 30 Myr and the 50% radius 
at r < 50 Myr. These features are due to the large number 
of BHs moving outward in the cluster on their way to escap- 
ing. The majority of these BHs have been removed from the 
cluster by ~ 50 Myr. After this point, the Lagrangian radii 
are still greatly inflated compared to those for Run 4. This 
is due to the extra kinetic energy received by the retained 
BHs in Run 5 at their creation. 

Fig. [26] shows the long-term evolution of the BH sub- 
system in Run 5. Because of the extra kinetic energy the 
retained BHs in this model receive at birth, they take sig- 
nificantly longer to sink to the cluster centre via dynamical 
friction than do the BHs in previous models. In addition, 
there are fewer BHs retained in Run 5, so once they have 
accumulated in the cluster core, they interact more infre- 
quently than in previous comparable models such as Run 4. 
The first BH binary does not form in Run 5 until r = 1350 
Myr, much later than in our previous models. Binary hard- 
ening, BH scattering and BH ejection subsequently begin; 
however, the rates of all these processes are considerably re- 
duced compared to previous simulations. The first ejection 
of a single BH after the formation of the first BH binary 
does not occur until r = 2050 Myr, while the first ejection 
of a BH binary does not occur until r = 4000 Myr. 

By the end of Run 5, at Tmax = 10.06 Gyr, only six BH 
binaries have been ejected. In common with earlier models, 
a significant population of BHs still remains in Run 5 at 
this point, consisting of 44 single BHs and one BH binary. 
The ejected BH binaries possess properties very similar to 
those observed for Run 4. Two have eccentricities in the 
range 0.8 < e < 0.95 while the remaining four have e < 
0.8. The closest ejected BH binary has separation ab — 7.55 
AU, while the least tightly bound has ai, = 61.0 AU. The 
mean mass of individual BHs in the ejected binaries is again 
greater than the ensemble average, at 10.65 M©. 

The reduced activity of the central BH subsystem in 
Run 5 compared with our other models explains the some- 
what slower expansion of the core radius in this simulation. 
Despite this, the evolution of rc/r^^i is very similar to earlier 
models with retained BH populations. Once the late, pro- 
longed phase of core expansion begins in Run 5 (i.e., after 
r ~ 1400 Myr), the locus traced by rc/rh,i lies exactly on top 
of that for Run 4, reaching rc/rh,i ~ 0.8 by the end of the 
simulation. This indicates that despite the reduced heating 
rate due to the BH population (and hence slower expansion 
of Tc), the distribution of this heating within the cluster is 
similar to that for Runs with larger numbers of BHs. Over- 
all, Run 5 clearly demonstrates that complete BH retention 
is not necessary for significant and prolonged core expan- 
sion to occur - even with half the number of retained BHs 
as Run 4, Run 5 still reaches the upper right-hand corner of 
the radius-age plane after ~ 10 Gyr of evolution. 
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4.5 Run 6: Can neutron stars replace BHs? 

Finally, we computed one further model, designed to inves- 
tigate whether the influence of a population of BHs is nec- 
essary for prolonged core expansion to develop in a cluster, 
or whether such expansion can also result from similar dy- 
namical processes involving larger numbers of less massive 
stellar remnants such as neutron stars (NSs). To this end 
we set up the new simulation, named Run 6, so that it was 
initially identical to Run 1 - that is, possessing no primor- 
dial mass segregation and retaining no BHs. However, unlike 
Run 1 where NSs were formed with large kicks so that the 
NS retention fraction /ns = 0, in Run 6 we set the natal NS 
kicks to be zero in order to achieve the extreme case that 
/ns = 1. In all, 425 NSs are formed in Run 6 from supernova 
explosions occurring between r « 10 — 43 Myr. The masses 
of these NSs lie in the range 1.30 5S mNS 5; 2.32 M0. We ex- 
tended the duration of Run 6 to be as long as that for Run 
1 (i.e., Tniax ~ 20 Gyr), to enable a comparison between the 
full long-term development of the two models. 

The evolution of Run 6 on the radius-age plane may 
be seen in Fig. 1271 Clearly, the retention of a large num- 
ber of NSs in this model does not result in core expansion 
at any stage during the lifetime of the cluster. In fact, the 
evolution is remarkably similar to that of Run 1, with the 
cluster undergoing many Gyr of gradual contraction before 
entering the core collapse phase. By the end of the early 
rapid mass-loss period at roughly r ~ 100 Myr, the median 
relaxation time in Run 6 is very similar to that in Run 1 
at the same age - i.e., trh ~ 2 Gyr. In the absence of any 
retained BHs, the NSs are the most massive objects in the 
system and hence begin to sink to the cluster centre on a 
time-scale of ~ (m*/mNs)ir-h ~ 500 Myr. However, the NSs 
are not very many times more massive than the otherwise 
most massive stars in the cluster, and so the central density 
of NSs never exceeds that of the other cluster members by 
a sufficient degree that the NS subsystem is unstable to a 
runaway collapse (|Spitzerlll987l . Eq. 3-55). Hence, the NS 
subsystem evolves quite differently to the BH subsystems 
in our previous models, which did become unstable to run- 
away collapse. No NS binaries are formed in the central core, 
and consequently, widespread scattering and ejection of NSs 
does not occur. As a result. Run 6 proceeds towards core 
collapse rather than undergoing prolonged core expansion. 

From Fig. [27] it is evident that Run 6 enters the core 
collapse phase at a significantly earlier time than does Run 
1 - observationally, the point of deepest collapse (smallest 
core radius) occurs aX t 'p^ 12.8 Gyr in Run 6, compared 
with T « 17.4 Gyr in Run 1. At any given age, the median 
relaxation time in Run 6 is very similar to that in Run 1, 
so that the point of deepest collapse in Run 6 occurs after 
significantly fewer integrated median relaxation times than 
in Run 1 - 4 AO trh as opposed to 8.37 trh- More enlightening 
is to examine the relaxation time in the central core of each 
cluster, tro oc Vsr^m~^ , where Vs is the velocity scale in 
the core and m,o is the mean mass of all t he particles i n 
thermal equilibrium in the central parts (e.g.. lMevlanlll987[ ). 
Calculating for each model the integrated number of central 
relaxation times which have elapsed by the time the point 
of deepest collapse occurs, the two values are within ~ 10 
per cent of each other. The central relaxation time in Run 
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Figure 27. Core radius evolution of A^-body Run 6, with the evo- 
lution of Run 1 plotted for comparison. These two models possess 
identical initial conditions — neither has any primordial mass seg- 
regation, and /bh = in both. The only difference between them 
is that in Run 6 neutron stars are formed with no natal kicks so 
that /ns = li whereas in Run 1 they are formed with large natal 
kicks so that /^g = 0. Hence, Run 6 retains 425 neutron stars, 
which are not present in Run 1. Unlike a BH population, the NS 
population in Run 6 does not lead to core expansion, but does 
cause the cluster to enter the core collapse phase at an earlier age. 



6 is generally much shorter than in Run 1, due to the larger 
mean mass of the centralmost objects (predominantly NSs) . 

It is also evident from Fig. 1271 that during collapse, the 
smallest observed core radius in Run 6 is larger than the 
smallest observed core radius in Run 1. This is due to the 
luminosity cut-offs inherent in the calculation of r^- In Run 
1, the stars contributing most of the light for the calculation 
of rc are also among the most massive remaining members, 
and hence typically reside in the innermost cluster regions 
during the core collapse phase. In Run 6 however, the cen- 
tralmost members are the NSs, which do not contribute light 
for the calculation of rc- The most luminous stars in Run 
6 therefore appear to constitute a more widely distributed 
"core" during the late phases of evolution than do those in 
Run 1. The larger core radius during collapse in Run 6 is 
also refiected in the evolution of rc/rh,i for this model. While 
the behaviour of this ratio is very similar to that for Run 1, 
during collapse rc/rh,i oscillates around ~ 0.2 rather than 
the smaller values {rc/rh,i < 0.05) observed for Run 1. 

Just as in Run 1, the point of deepest collapse in Run 
6 coincides with a spate of binary formation in the core. 
This time, the binaries generally possess at least one NS 
member; several of them are NS-NS binaries. These are the 
first binary objects involving NSs to be created in Run 6 - 
no such objects are formed at earlier times in this model. 
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Figure 28. Our full suite of eight A^-body Runs plotted together for comparative purposes. In combination, the two core expansion 
processes described in this paper lead to a wide variety of evolutionary paths on the versus age plane, which fully cover the observed 
distribution of massive Magellanic Cloud star clusters. 



5 DISCUSSION AND SUMMARY 

In this paper we have presented an ensemble of eight large- 
scale A''-body simulations aimed at directly modelling the 
evolution of massive star clusters like those observed in the 
Magellanic Clouds. Figure [2S] shows the core radius evolu- 
tion of all eight models plotted on the same set of axes, 
for the purposes of direct comparison. Using these models 
we have identified two physical processes which can occur 
in such clusters and result in substantial core expansion - 
mass-loss due to rapid stellar evolution in a cluster which 
is mass-segregated or otherwise centrally concentrated, and 
heating due to a significant population of retained stellar- 
mass (~ 10 M0) BHs formed in the supernova explosions of 
the most massive cluster stars. These two processes operate 
over different time-scales and at different stages in a clus- 
ter's life. The former only occurs during the first ~ 100 — 200 
Myr after the formation of a cluster, when massive stars are 
still present. These evolve rapidly, losing a large fraction of 
their mass as they do so. The latter begins, at the earliest, 
several hundred Myr after the formation of the cluster but 
may remain active for at least a Hubble time beyond this 
starting point. In combination, these two processes can lead 
to a wide variety of evolutionary paths on the core-radius 
versus age plane, which fully cover the observed distribution 
of massive star clusters. They therefore define a physically- 



motivated dynamical explanation for the radius-age trend 
observed for the star cluster systems belonging to the Mag- 
ellanic Clouds. 

Our A'^-body modelling has allowed us to examine in 
more detail the behaviour of each of these core-expansion 
mechanisms. As stated above, the phase of severe mass- 
loss due to rapid stellar evolution is mostly complete by 
~ 100 — 200 Myr into a cluster's life, by which time all the 
most massive stars in the cluster have expired. Mass-loss due 
to stellar evolution still occurs after this point; however, it 
is from much less massive stars and therefore proceeds at a 
far more gradual rate without noticeably affecting the core 
radius evolution of the host cluster. Our simulations show 
that the amount of observed core expansion in a cluster due 
to the early mass-loss phase depends on both the degree of 
primordial mass segregation in the cluster, and the amount 
of mass lost in relation to the total cluster mass. In models 
where the former parameter is held constant and the latter 
parameter is varied (e.g., Runs 3, 4, and 5), the cluster los- 
ing the most mass expands the fastest. In models where the 
latter parameter is held constant and the former parame- 
ter is varied (e.g.. Runs 2, 4, 4a, and 4b), the cluster with 
the most significant degree of primordial mass segregation 
expands the fastest. Furthermore, the early rapid phase of 
mass-loss does not cause any significant core expansion in 
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our models unless the mass-loss is centrally concentrated 
- models which do not possess any primordial mass seg- 
regation exhibit essentially no early expansion. In models 
where early expansion occurs, the ratio of the core radius 
to half-light radius rc/rh,i increases significantly. This is in 
contrast to models which do not undergo early expansion, 
where r^jru^i remains fairly constant with time. Inflated val- 
ues of rc/rh,i may remain for > 10 Gyr in some clusters (cf. 
Run 3), since the central and median relaxation times in 
these objects become rather long. 

Since the amount of mass lost in the early rapid stellar 
evolution phase is an important contributor to the amount 
a cluster core expands during this phase, the expansion is 
effectively regulated by the cluster's stellar IMF, modulated 
by second order effects such as BH retention. A steep IMF re- 
sults in proportionally few high-mass members and hence a 
small amount of early expansion, whereas a flat IMF results 
in proportionally many high-mass members and hence a 
large amount of early expansion. In principle, therefore, sig- 
nificant inter-cluster IMF variations could lead to a variety 
of dramatically different paths across the radius-age plane 
at ages t < 200 Myr, and consequently in duce a large sprea d 
in the observed cluster distribution (e.g.. lElson et al.lll989h . 
However, there is an increasing body of evidence that large- 
scale variations in the stellar I MF between Magellanic Cloud 
cluste rs are not present (e.g., iKroupal I2001I : Ide Griis et al.l 
l2002d ). This in turn suggests that inter-cluster variations in 
the degree of primordial mass segregation or central concen- 
tration may be the primary driver of the observed spread in 
the radius-age distribution at young ages. The sharp upper 
envelope of the distribution at ages less than a few hun- 
dred Myr therefore most likely reflects an upper limit to 
the degree of primordial mass segregation present in Magel- 
lanic Cloud clusters. Indeed, our model with an IMF match- 
ing that gene rally observed for Magellanic Cloud clusters 
(|Kroupall200ll \ and an initial structure (including mass seg- 
regation) matching that observed for R136 (which is the 
most extreme young object presently observed in the Mag- 
ellanic Clouds) evolves along a path closely matching the 
upper envelope of the radius-age distribution at early times. 

One important process which can affect the early evo- 
lution of a massive star cluster, but which is not included 
in our A'^-body models, is the expulsion of gas left over 
from the star formation process. This residual gas is re- 
moved from the cluster within the flrst ~ 10 Myr due to 
the combined effect of massive stellar winds and supernova 
explosions. Just as with the early mass-loss due to stellar 
evolution, mass-loss due to gas expulsion can cause clus- 
ter core expansion, although ty pically on an even shorter 
time-scale of ^ 10 - 20 My r fe.g.. lBastian fc Goodwinll2006l : 
ICoodwin fc BastianI 120061 ). The larger the mass of expelled 
gas, the larger the degree of core expansion; if the mass of 
expelled gas is too great, the cluster may become unbound. 
In cases where the star formation efficiency is relatively high, 
and in the absence of sustained mass-loss due to stellar evo- 
lution, clusters soon settle into new equilibrium states with 
core radii generally not much larger than their initial values 
(see lCoodwin fc Ba stian 2006, Fig. 2). Therefore, gas expul- 
sion may be affecting the core radius evolution of clusters 
younger than ~ 50 Myr in our radius-age plot, but is unlikely 
to be of relevance to cluster evolution on longer time-scales. 
To correctly model the effects of gas loss occurring in com- 



bination with early stellar mass loss on the evolution of the 
various types of clusters studied in the present work will re- 
quire more sophisticated codes than are presently available. 

Core expansion due to the dynamical influence of a pop- 
ulation of retained stellar-mass BHs in a cluster occurs over 
a much longer timescale than that due to early mass loss. 
Our simulations show that the BH population in a cluster 
only induces core expansion once the BHs have accumulated 
in a sufficiently dense central subsystem that BH binaries are 
created. These binaries are the catalyst for core expansion, 
since it is the interactions between them and other single 
and binary BHs which lead to BH scattering and ejection, 
and subsequent heating of the central cluster regions. We do 
not observe core expansion due to the BH population prior 
to the formation of BH binaries in any of our simulations. 
In particular, the sinking and accumulation of BHs in the 
core does not appear to affect the observed core radius. 

In our models, the time at which the first BH binaries 
are formed is relatively independent of the early evolution 
of the cluster. Models which are identical but for widely 
varying degrees of primordial mass segregation and hence 
widely varying amounts of early expansion (Runs 2, 4, 4a, 
and 4b) all form their first BH binaries at ages of roughly 
~ 500—600 Myr. Even though BHs are preferentially formed 
very centrally in a model with a significant degree of pri- 
mordial mass segregation compared to a model with no pri- 
mordial mass segregation, the former object undergoes very 
significant early expansion compared to the latter. The BH 
subsystem does not escape this expansion, and by r ~ 200 
Myr it has roughly the same distribution within the cluster 
as does the BH subsystem in the initially non-segregated 
model. The subsequent evolution of the two BH populations 
is very similar. The time of formation of the first BH bina- 
ries is, however, strongly sensitive to the natal kicks received 
by the BHs at formation. In the case of non-zero kicks, re- 
tained BHs take longer to accumulate in the cluster centre 
than in the case of no kicks, due to the extra kinetic en- 
ergy they receive at birth. In addition, non-zero natal kicks 
generally result in the expulsion of some fraction of the BH 
population, leading to a smaller retained BH subsystem and 
a smaller probability per unit time of BH binary formation. 
Although our modelling did not test it, the time of formation 
of the first BH binaries is also expected to be sensitive to the 
mean BH mass. More massive BHs will sink to the cluster 
centre much more rapidly than less massive BHs, and hence 
form a dense central core at a significantly earlier time. 

Once BH binaries have formed in a cluster and core ex- 
pansion begins, the rate of expansion is dependent on the 
number of BHs in the cluster. A cluster with fewer BHs ex- 
pands more slowly than an otherwise identical cluster with 
more BHs (cf. Runs 4 and 5) . This is because the interaction 
rate between BHs in the cluster centre is much lower for the 
model with the smaller number of BHs, so that fewer BHs 
are scattered and ejected per unit time and the rate of heat- 
ing of the cluster is reduced. The interaction rate between 
BHs in the cluster centre is also apparently sensitive to the 
density of the surrounding stellar core - it is significantly 
reduced in lower density clusters (cf. Runs 2 and 4). These 
observations have important implications for the survivabil- 
ity of BH subsystems within clusters. As the number of BHs 
in a cluster decreases due to the ejection of BHs after close 
encounters, and the central density of the cluster decreases 
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due to the expansion of its core, the interaction rate between 
BHs in the cluster centre also decreases. This in turn results 
in a lower BH ejection rate, allowing the BH population in a 
cluster to survive much longer than previously believed. All 
our long-duration simulations with retained BHs still possess 
a sizeable BH population after a Hubble time of evolution. 
As a result, some degree of core expansion is still occurring 
in these models at late times. 

We emphasize that even though most of our models 
examine the scenario where all BHs are retained in a cluster, 
such an extreme case is not required for core expansion to 
occur. We still observe significant expansion in the more 
moderate case of ~ 50 per cent retention, although the rate 
of expansion is reduced due to the factors outlined above. 

It is also worth emphasising that while rapid mass-loss 
due to stellar evolution is the dominant cluster core expan- 
sion process at early times (r < 200 Myr), that expansion 
ceases as the mass-loss rate slows. This process therefore 
cannot drive the full observed radius-age distribution, which 
still exhibits a significantly increasing spread in core size at 
much later times. A cluster which has expanded during its 
first few hundred Myr of evolution, but which has not re- 
tained a BH population sufficient to induce additional late- 
time expansion, begins to contract again as two-body relax- 
ation processes take over (Run 3). Our models only achieve 
core sizes matching those observed for the most extended 
~ 10 Gyr old Magellanic Cloud clusters if expansion due 
to a retained BH population also occurs. This long-term 
expansion cannot be reproduced by other types of stellar 
remnants such as neutron stars (Run 6). Such remnants are 
not of high enough mass to accumulate in a central subsys- 
tem of sufficient density to allow frequent formation of the 
binary objects which drive the cluster heating. 

The ratio rc/rh,i evolves very similarly in all of our mod- 
els where core expansion due to a retained population of BHs 
occurs, once this phase has started. By r ~ 10 Gyr, the ra- 
tio approaches a large value of rc/rh,i ~ 0.8, comparable to 
the largest values observed for old Magellanic Cloud clus- 
ters, and Galactic globular clusters. This is irrespective of 
the early evolution of a cluster (i.e., whether expansion due 
to early mass-loss occurs or not), the time of onset of the 
expansion due to BHs, and the subsequent rate of this ex- 
pansion. As described in Section |4j these o bservations ar e 
compatible with those presented recently bv lHurlevI l|2007h . 

Several other mechanisms are known to be able to sus- 
tain large or expanding cores in massive clusters. For exam- 
ple, the presence of primordial binary stars may stall core 
collapse, while the presence of a central IMB H may result in 
cluster expansion l|Baumgardt et all l2004al lbh . The heating 
effect of stellar-mass BHs, as considered in this paper, is far 
more efficient than the heating effect due to primordial bina- 
ries in comparable clusters. To transfer binding energy from 
primordial binaries to other cluster members requires fre- 
quent interactions and hence a dense environment. For most 
of their lives, Magellanic Cloud clusters are not sufficiently 
de nse, as has be en demonstrated from A*'-body modelling 
bv lMackevI (|2003l 1^. Furthermore, heating due to primordial 
binaries is self-regulated: a dense core will expand, reducing 
the interaction rate and switching the heating off until the 

* This Ph.D. Thesis can be supplied by ADM on request. 



core contracts again. Primordial binaries therefore cannot 
sustain the type of core expansion observed in our A'^-body 
models. It is more difficult to estimate the relative heating 
efficien cy of stellar mass BHs com pared to that of a central 
IMBH. iBaumgardt et al.l (|2004al lbl) display the evolution of 
the Lagrangian radii of their large A-body models, which do 
show significant expansion. However, it is difficult to disen- 
tangle the amount of heating due to mass-loss from stellar 
evolution from the amount due to the central I MBH. Based 
on the material presented bv lBaumgardt et al.l (|2004al lbh ■ we 
estimate that heating due to stellar-mass BHs is probably 
at least as efficient as that due to a central IMBH. 

The scenario outlined above as a dynamical explanation 
for the radius-age trend observed in the Magellanic Clouds 
requires significant variations in BH population size (i.e., in 
the BH retention fraction) between otherwise very similar 
clusters. Clusters which have developed very large core radii 
by the time they are > 12 Gyr old (e.g., the LMC clusters 
NGC 1841 and Reticulum) must have managed to retain a 
significant BH population. Conversely, clusters which have 
entered core collapse at late times (e.g., the LMC clusters 
NGC 2005 and 201 9) are unlikely to have retained very many 
BHs - for example. iHurlevI (|2007l ') showed that even one BH 
binary in a cluster can prevent the collapse of its core. 

There are a number of possibilities which could lead 
to inter-cluster variability in the BH population size. First, 
we note that the number of BH-forming stars in a given 
cluster is only a tiny fraction of the total number of stars in 
the cluster, so there will always be sampling-noise variations 
between clusters. In addition, the formation of a BH in a su- 
pernova explosion is sensitive to many features of the prior 
evolution of the progenitor star, in particular how much 
mass it loses as it evolves. Factors which introduce mass- 
loss variations, such as binarity, chemical inhomogeneities 
or a dispersion in stellar rotation, are therefore likely to fur- 
ther accentuate the stochastic fluctuations in BH population 
size between clusters. In principle, inter-cluster variations in 
the stellar IMF would also strongly affect BH population 
sizes; however, as we noted earlier, such variations are not 
observed. Observations do suggest that the maximum stellar 
mass in a cluster correlat es with the total cluster mass (e.g., 
IWeidner fc Kroupalliooil '). Hence, even if the steUar IMF is 
universal between clusters, smaller clusters will have a lower 
maximum stellar mass and thus fewer BHs relative to the 
total cluster mass than will larger clusters. 

Natal BH kicks are also critical. At present, these are 
poorly constrained both by theory and observation. Typical 
estimates lie in the range < iikick ^ 200 km s~^, with kicks 
of a f ew tens of kms~^ possibly favoured (e.g.. lWillems et al.l 
|2005| . and references therein). If BH natal kicks are indeed 
typically a few tens of kms~^ in magnitude, then they are 
roughly comparable to the escape velocity of a massive stel- 
lar cluster. In this case, the structure of the host cluster when 
the BHs are formed (i.e., before r « 10 Myr) can strongly 
affect the retention fraction. For example, BHs formed in a 
dense, strongly mass segregated cluster are more easily re- 
tained that BHs formed in a comparably massive but less 
dense, non-segregated cluster (e.g.. Fig. [23]). The retention 
fraction will also be sensitive to the overall initial mass of the 
cluster, as well as to factors which affect the very early evo- 
lution of the cluster such as residual gas expulsion. Stellar 
binarity may also play an important role. 
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It is interesting to note that theoretical models sug- 
gest BH formation to be a strongly sensitive function of 
metal abundance, in that BH production is apparently more 
frequent, and rriBH is typically greater for metal poor pro- 
genitor stars than for metal rich progenitor stars (see e.g., 
IZhang et al"]l2007l l. Hence, the BH populations formed in 
clusters of very different metallicities are likely to be quite 
distinct. The strong age-metallicity relationships observed 
in both the LMC and SMC mean that this factor prob- 
ably cannot cause differences between the BH population 
sizes in compact and extended LMC or SMC clusters of a 
given age, since such objects will have approximately equal 
metallicities. However, the LMC and SMC age-metallicity 
relationships do suggest that any BH populations forming 
in present-day young Magellanic Cloud clusters are likely to 
be quite different to those which may have existed in Mag- 
ellanic Cloud clusters that are now > 12 Gyr old. 

The possibility of large-scale and prolonged core expan- 
sion has important implications for the study of all massive 
star clusters, including the globular clusters in our Galaxy 
and others. Many such objects are at least an order of magni- 
tude more massive than the models presented in this paper. 
Even so, we expect the physical processes which we have 
described will still operate in larger systems. 

Neglecting any stochastic fluctuations between clusters, 
the number of BHs formed in a cluster is, to first order, de- 
pendent only on the stellar IMF and the minimum progen- 
itor mass. We do not expect these to change with cluster 
mass, so with all other parameters held constant, the mass 
fraction of BHs should remain the same for clusters of in- 
creasing mass. Similarly, the mass fraction lost due to rapid 
stellar evolution early on in a cluster's life should also re- 
main the same for clusters of increasing mass. Assuming 
that natal BH kicks are also not a function of cluster mass, 
the BH retention fraction should increase with increasing 
cluster mass, since it is more difficult to eject BHs from the 
deeper gravitational well of a more massive cluster. Overall, 
we therefore expect the relative size of retained BH popula- 
tions should be larger for more massive clusters. Given the 
above, more massive clusters have at least the same poten- 
tial for core expansion as do less massive clusters. 

In terms of the dynamics of the expansion, the central 
densities of the model clusters we have studied in this pa- 
per are directly comparable to the central densities of more 
massive objects such as globular clusters. The central and 
median relaxation times in our models are also commen- 
surate with those calculated for typical globular clusters. 
Given this, we expect similar dynamical processes to oper- 
ate on similar time-scales in clusters larger than our present 
models, so that early mass-loss and BH heating will both 
still be effective at inducing core expansion. The main dif- 
ference is that it becomes more difficult to eject BHs as the 
cluster mass increases. Therefore, the mean time a BH re- 
mains in a cluster will increase with the total mass. This will 
increase the potential of each BH to heat the cluster through 
additional scattering-sinking cycles, and will allow BH bina- 
ries to harden further than they would do in a less massive 
object. BH heating in more massive clusters is hence likely 
to be even more efficient than it is in less massive clusters. 

We therefore predict that some degree of core expansion 
is possible in any massive stellar cluster due to the processes 
outlined in this paper, irrespective of the total mass of the 



cluster. In many aspects of star cluster research, this possi- 
bility is not usually considered. However, under a wide vari- 
ety of circumstances, it could have an important effect on the 
problem under consideration. As a simple example, it is well 
known that diffuse extended clusters are far more suscepti- 
ble to tidal disruption than are compact clusters. Prolonged 
core expansion in clusters could result in many more such 
diffuse extended objects in a given population than would 
otherwise be expected. This possibility is vital to incorporate 
into modelling where destruction rates are important, such 
as studies designed to investig ate the evolution of the globu- 
lar cluster mass function (e.g.. lFall fc Zhan3l200ll ). the past 
and future d issolution of globular clu sters in the Galactic 
system (e.g., iGnedin fc Ostrikeij|l997i ). or whether the su- 
per star clusters observed in starburst galaxies will eventu- 
ally evolve into objects resembling classical globular clusters 
(e.g., de Griis & Parmcnticr 2007). 

Another example involves the measurement of dynam- 
ical mass estimates for young massive clusters in external 
galaxies. Such measurement s are sometimes used to infe r the 
stellar IMF in such clusters. iBastian fc GoodwinI (|2006l ) and 
iGoodwin fc BastianI (j2006l ) have demonstrated that very 
young clusters (r < 50 Myr) may be out of virial equilibrium 
due to the rapid expulsion of residual gas, so that dynamical 
mass measurements assuming virial equilibrium may be in 
error. Our modelling has shown that significant core expan- 
sion due to stellar evolution occurs on a timescale close to 
~ 100 Myr. Researchers should be aware of this additional 
possibility when evaluating the properties of young clusters 
in external galaxies, although we note that it is not yet clear 
to what extent any signal due to such expansion would man- 
ifest in integrated cluster spectra. This is an avenue worthy 
of further investigation. 

As a final example, consider the cluster half-mass ra- 
dius, rh- This quantity is often assumed to be relatively 
stable for much of a cluster's life (cf. Fig. O, and is hence 
some times used to infer information a bout cluster formation 
(e.g.. Ivan den Bergh fc Mackevll2004l '). However, if a cluster 
undergoes prolonged core expansion, is certainly not a 
stable quantity (Fig. [8} . Caution should therefore be exer- 
cised in the use of this parameter. 

The possibility of core expansion may also help explain 
the properties of some of the more exotic star clusters dis- 
covered in recent years - for example, t he "faint fuzzies" 
unco vered in several lenticular galaxies (|Brodie fc LarsenI 
I2002D the luminous extended clusters found in the halo of 
M31 (|Huxor et al.ll2005l : lMa"ckev et al.ll2006l ). and the diffuse 
star clusters located i n a number of Virgo early- type galax- 
ies (|Peng et al.l |2006| ). All these newly-discovered clusters 
possess unusually extended structures compared to those of 
classical globular clusters. Core expansion, particularly the 
prolonged variety due to retained BHs, may offer a viable 
formation channel for these objects. 

We conclude with a note on the possibility of test- 
ing observationally our prediction of retained populations 
of stellar-mass BHs in some massive star clusters. While 
these BHs cannot be observed directly, there are a num- 
ber of means by which their presence might be inferred in 
a cluster. One possibility is that a BH in a close binary 
with an evolved star is likely to be an X-ray source, as 
the star overffows its Roche limit and transfers gas onto 
the BH. Such BH X-ray sources are known in the field (see 
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e.g., ICasared |2006|) and one is known in an extra-Galacti c 
globular cluster l|Maccarone et al] l2007l : HeEfeFai] |20o3) ; 

however, none have been found in Galactic globular clusters 
IVerbunt fc Lewi n 20061 . From our modelling, we know that 
clusters which do retain significant BH populations are, for 
most of their lives, objects in which the timescale for en- 
counters between BHs and stars is very long (due to the low 
stellar density in the extended core), but the timescale for 
encounters between BHs is comparatively short. Ifence the 
creation of long-lived BH-star binaries is rare - we did not 
observe any such objects in our simulations. It is therefore 
unsurprising that no BH X-ray sources are known in the 
Galactic globular cluster system, and only one is known in 
an external cluster. 

The most promising means of inferring the presence of 
a BH population in a cluster is through the dynamical ef- 
fect it causes on the stellar component of the cluster. As we 
demonstrated in Section r4.1l unlike in the case of an IMBH, 
a significant stellar density and velocity cusp does not de- 
velop about the compact central BH subsystem. None the 
less, the effect of the central concentration of unseen mass 
should be evident in the stellar motions - the velocity disper- 
sion of the cluster should be larger than is to be expected 
solely from the observed luminous mass. Observations to 
test this will be difficult, primarily because the target clus- 
ters should be extended, diffuse objects with low velocity 
dispersions. In addition, many will have relatively few tar- 
gets suitable for spectroscopic radial velocity measurements, 
such as luminous red giants. Even so, it may be possible to 
make sufficiently precise observations with presently avail- 
able 8 — 10 m-class facilities. 

Finally, it seems likely that at least some BH bina- 
ries ejected from very massive clusters will merge on a 
timescale < 12 Gyr, and will therefore be sources of grav- 
itational radiation detectable by interferometers such as 
LIGO, and in future, LISA. This possibility has previously 
been investigated in more detail b y other authors (see e.g., 
iPortegies Zwart fc McMillanll2000l '). As described in Sections 
I4.1l and l4.2l the BH binaries ejected from our model clusters 
will not merge on a timescale < 12 Gyr; however several 
have orbital parameters not far from the required thresh- 
old. A subset of the BH binaries ejected from more massive 
clusters than those studied here would almost certainly have 
orbital parameters well within this threshold. 
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APPENDIX A: ANALYTIC PROPERTIES OF THE EPF FAMILY OF MODELS 

In this Appendix we present analytic expressions for the properties of a number of members of the general family of EPF 
models. As demonstrated in the present work, with the recent rapid increase in computing power and software sophistication, 
and hence the size and degree of realism feasible for A''-body simulations, it may occur that such models represent more 
appropriate initial conditions for a given scenario than do the frequently adopted Plummer spheres or King models. With 
the formulae presented below, it is reasonably straightforward to develop procedures such as that described in Section [2. II to 
generate the desired initial conditions. 



Al General properties 

The EFF models, after lElson. Fall fc Freeinaiil (|l987l ). are a family of models originally defined through empirical fitting to 
the observed surface brightness profiles of young massive star clusters in the LMC. These objects do not exhibit tidal cut-offs 
in their outer regions, and are therefore most appropriately described by projected three-parameter models of the form: 



' 1 



^l{r^) = ^0 (^1 + -I j , (Al) 

where is the projected radius, is the central surface brightness, a is the scale radius, and 7 represents the power-law 
fall-off of the profile at large radii. These models can easily be deprojected to obtain the three-dimensional density: 

pW = Po 1 + - where Po ^ ^ w% v (A2) 

In the above equation, F is a standard gamma function. Since /io is the central surface brightness, here po is the central 
luminosity density - to obtain the central mass density it is necessary to multiply by the global mass-to-light ratio T. It can 
be seen that the three-dimensional density has exactly the same functional form as the projected densi ty, but with index 
7 -I- 1. By comparison with the more general spherically symmetric family of (a, j3, 5) models described by [Zhaol il996 |lPI, it is 
straightforward to see that the EFF models represent the subset with (a, /3, S) = (^,7-1- 1, 0) and break radi us r — a. 

Assuming now that po is a mass density, substituting Ea. lA2l into Eq. 2-22 from lBinnev fc Tremainj (|l987l l and integrating 
yields a general expression for the gravitational potential of EFF models: 



4 

>!> = --TvGpo 



(t-1) 

3aW, ^ 2 r. 7 + 1 5 



(A3) 



where 2-Fi(a, b; c; z) is Gauss's hypergeometric function. 

Similarly, the enclosed mass (or luminosity if po is a luminosity density) as a function of radius can be derived by 
integrating Eg. IA2I 



M(r) = 4n / p(r')r'^dr' = ^npor' ^J'W | l±i; | - ^ ) . (A4) 







In the limit where r —+ 00, this expression converges only if 7 > 2. The asymptotic mass is given by Moo ~ /{'y — 2). 

Finally, rearranging and int egrating the Jeans equations for a steady-state, spherically symmetric, non-rotating cluster 
(i.e., iBinnev fc Tremaind (|l987l ) Eq. 4-54) yields a general expression for the radial dependence of the isotropic velocity 
dispersion: 

which can, in principle, be evaluated numerically for all 7. However, for integer values of 7, the hypergeometric functions in 
Eos. lASl and IA4l reduce to elementary functions, resulting in straightforward analytic expressions for M(r) and o"^(r) which 
may be wri tten into t h e com puter code to generate initial conditions. The best known of the analytic su bset is the case w hen 
7 = 4- the IPlummeij (|l91lh sphere. The properties of this m odel have been investig ated extensively bv lDeionghel (|l987l ) for 
mass-foUows-light scenarios, while the more general study of I Wilkinson et al] l|2002l ) includes the possibility of a dark halo. 
Below, we consider the less-well studied cases of 7 = 3, 5, and 6. 



A2 The 7 = 3 case 

The case where 7 = 3 has been studied in passing as special cases of a general ellipsoidal form bv lde Zeeuwl (|l985al lbh. who 
labelled the profiles "perfect spheres". This particular case is of interest since it represents the EFF model with analytic 
expressions for M(r) and (T^(r-) which, in projection, is closest in form to the observed profiles of young massive star clusters 



^ IZhaol lll996l ) labelled these {a, 0, 7) models; however, we alter his nomenclature here to avoid ambiguity in the definition of 7. 
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in the LMC and SMC. lElson et al.l l|l987^ found a median value of 7 = 2.6 and a range 2.2 < 7 < 3.2 for tiieir ten young LMC 
clusters, while the 18 young LMC and SMC clusters plotted in Fig.[3]in the present paper cover the range 2.05 ^ 7 ^ 3.79 and 
have a median value 7 = 2.67. We started all the A'^-body simulations described in the present work with initial conditions 
generated from the following equations. 

When 7 = 3, the Gauss hypergeometric function in Egs. IA3l and IA4| reduces to: 
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Substituting into Eq. IA4I the enclosed mass M{r) is then given by: 



M{r) — 2npo a"^ ^arctan 



1 + ^ 



while carrying out the integration in Eg. I A5I provides the isotropic velocity dispersion 
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4 ( ar + {a^ + r^) arctan 



4a'' + 3ar^ + 3r(a^ + r^) arctan 



(A6) 



(A7) 



(A8) 



r{a^ + 7-2)2 
A3 Steeper cases: 7 = 5 and 7 = 6 

Cluster models with steep density fall-offs do not seem to have been well studied in the literature, if at all. For this reason, 
the properties of two such models are derived here. These 7 = 5 and 7 = 6 cases do co rrespond to real objects - o ld globular 
clusters can have observed brightness profiles which fall off this steeply. For example, iMackev fc Gilmord (|2003al lblld) found 
the LMC globular cluster NGC 1841 to have 7 = 4.55 ± 0.61, the SMC cluster NGC 339 to have 7 = 5.21 ± 0.99, and the 
diffuse cluster 1 in the Fornax dSph to have 7 = 7.52 ± 0.64. ft may well be desirable in future to model the evolution of 
clusters such as these. In this case, the equations below will allow suitable initial conditions to be simply constructed. 
If 7 = 5, the hypergeometric function in Eqs. IA3I and IA4I reduces to: 
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Substituting into Eq. IA4I as before, yields the enclosed mass M{r): 
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while evaluating Eq. IA5I provides the isotropic velocity dispersion as a function of radius: 
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The 7 = 6 case is of particular interest as its properties are comparable in simplicity to those of the widely used Plummer 
sphere (7 = 4). With 7 = 6, the hypergeometric function in Egs. I A3l and I A4l can be written: 

(A12) 

(A13) 
(A14) 



2F1 

which leads to the following straightforward expressions for the enclosed mass and isotropic velocity dispersion: 
M(r) 
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In principle, it is possible to continue deriving similar analytic expressions for increasing integer values of 7; however we 
note that the expressions become correspondingly more complicated as 7 increases. 



